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Abstract 

It  was  believed  until  very  recently  that  a  near-equatorial  satellite  would  always  keep 
up  with  the  planet’s  equator  (with  oscillations  in  inclination,  but  without  a  secular 
drift).  As  explained  in  Efroimsky  and  Goldreich  (2004),  this  misconception  originated 
from  a  wrong  interpretation  of  a  (mathematically  correct)  result  obtained  in  terms 
of  non-osculating  orbital  elements.  A  similar  analysis  carried  out  in  the  language  of 
osculating  elements  will  endow  the  planetary  equations  with  some  extra  terms  caused 
by  the  planet’s  obliquity  change.  Some  of  these  terms  will  be  nontrivial,  in  that 
they  will  not  be  amendments  to  the  disturbing  function.  Due  to  the  extra  terms, 
the  variations  of  a  planet’s  obliquity  may  cause  a  secular  drift  of  its  satellite  orbit 
inclination.  In  this  article  we  set  out  the  analytical  formalism  for  our  study  of  this  drift. 
We  demonstrate  that,  in  the  case  of  uniform  precession,  the  drift  will  be  extremely  slow, 
because  the  first-order  terms  responsible  for  the  drift  will  be  short-period  and,  thus, 
will  have  vanishing  orbital  averages  (as  anticipated  40  years  ago  by  Peter  Goldreich), 
while  the  secular  terms  will  be  of  the  second  order  only.  However,  it  turns  out  that 
variations  of  the  planetary  precession  make  the  first-order  terms  secular.  For  example, 
the  planetary  nutations  will  resonate  with  the  satellite’s  orbital  frequency  and,  thereby, 
may  instigate  a  secular  drift.  A  detailed  study  of  this  process  will  be  offered  in  the 
subsequent  publication,  while  here  we  work  out  the  required  mathematical  formalism 
and  point  out  the  key  aspects  of  the  dynamics. 
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1  Physical  motivation  and  the  statement  of  purpose 

Ward  (1973,  1974)  noted  that  the  obliquity  of  Mars  may  have  suffered  large-angle  motions 
at  long  time  scales.  Later,  Laskar  and  Robutel  (1993)  and  Touma  and  Wisdom  (1994) 
demonstrated  that  these  motions  may  have  been  chaotic.  This  would  cause  severe  climate 
variations  and  have  major  consequences  for  development  of  life. 

It  is  a  customary  assumption  that  a  near-equatorial  satellite  of  an  oblate  planet  would 
always  keep  up  with  the  planet’s  equator  (with  only  small  oscillations  of  the  orbit  inclina¬ 
tion)  provided  the  obliquity  changes  are  sufficiently  slow  (Goldreich  1965,  Kinoshita  1993). 
As  demonstrated  in  Efroimsky  and  Goldreich  (2004),  this  belief  stems  from  a  calculation 
performed  in  the  language  of  non-osculating  orbital  elements.  A  similar  analysis  carried 
out  in  terms  of  osculating  elements  will  contain  hitherto  overlooked  extra  terms  entailed  by 
the  planet’s  obliquity  variations.  These  terms  (emerging  already  in  the  first  order  over  the 
precession-caused  perturbation)  will  cause  a  secular  angular  drift  of  the  satellite  orbit  away 
from  the  planetary  equator. 

The  existence  of  Phobos  and  Deimos,  and  the  ability  of  Mars  to  keep  them  close  to 
its  equatorial  plane  during  obliquity  variations  sets  constraints  on  the  obliquity  variation 
amplitude  and  rate.  Our  eventual  goal  is  to  establish  such  constraints.  If  the  satellites’ 
secular  inclination  drifts  are  slow  enough  that  the  satellites  stay  close  to  Mars’  equator  during 
its  obliquity  changes  through  billions  of  years,  then  the  rigid-planet  non-dissipative  models 
used  by  Ward,  Laskar  &  Robutel,  and  Touma  &  Wisdom  will  get  a  totally  independent 
confirmation.  If  the  obliquity-change-caused  inclination  drifts  are  too  fast  (fast  enough  that 
within  a  billion  or  several  billions  of  years  the  satellites  get  driven  away  from  Mars’  equatorial 
plane),  then  the  inelastic  dissipation  and  planetary  structure  must  play  a  larger  role  than 
previously  assumed. 

Having  this  big  motivation  in  mind,  we  restrict  the  current  article  to  building  the  required 
mathematical  background:  we  study  the  obliquity-variation-caused  terms  in  the  planetary 
equations,  calculate  their  secular  components  and  point  out  the  resonant  coupling  emerg¬ 
ing  between  the  satellite’s  orbiting  frequency  and  certain  frequencies  in  the  planet  axis’ 
precession.  A  more  thorough  investigation  of  this  interaction  will  be  left  for  our  next  paper. 

2  Mathematical  preliminaries: 
osculating  elements  vs  orbital  elements 

Whenever  one  embarks  on  integrating  a  satellite  orbit  and  wants  to  take  into  account  direc¬ 
tion  variations  of  the  planet’s  spin,  it  is  most  natural  to  carry  out  this  work  in  a  co-precessing 
coordinate  system.  This  always  yields  orbital  elements  defined  in  the  said  frame  and,  there¬ 
fore,  ready  for  immediate  physical  interpretation  by  a  planet-based  observer.  A  well  camou¬ 
flaged  pitfall  of  this  approach  lies  in  that  these  orbital  elements  may  come  out  non-osculating, 
i.e. ,  that  instantaneous  ellipses  (or  hyperbolae)  parametrised  by  these  elements  will  not  be 
tangent  to  the  physical  trajectory  as  seen  in  the  said  frame  of  reference. 

2.1  The  osculation  condition  and  alternatives  to  it 

An  instantaneous  orbit  is  parametrised  with  six  independent  Keplerian  parameters.  These 
include  the  three  Eulerian  angles  i,  u>,  Q  which  define  the  orientation  of  the  instantaneous 
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orbital  plane,  the  eccentricity  and  semimajor  axis  of  the  instantaneous  ellipse  or  hyperbola, 
and  the  mean  anomaly  at  an  epoch,  M0  ,  which  determines  the  initial  position  of  the  body. 
Often  are  employed  sets  of  other  variables.  The  variables  are  always  six  in  number  and  are 
some  functions  of  the  Keplerian  elements.  These  are,  for  example,  the  Delaunay  set,  the 
Jacobi  set,  and  two  sets  offered  by  Poincare.  More  exotic  is  the  Hill  set  which  includes  the 
true  anomaly  as  one  of  the  elements  (Hill  1913). 

Systems  of  planetary  equations  in  terms  of  the  above  sets  of  parameters  may  be  de¬ 
rived  not  only  through  the  variation-of-parameters  (VOP)  method  but  also  by  means  of  the 
Hamilton-Jacobi  one.  However,  the  latter  approach,  being  fine  and  elegant,  lacks  the  power 
instilled  into  the  direct  VOP  technique  and  cannot  account  for  the  gauge  invariance  of  the 
N-body  problem  (Efroimsky  2002a, b),  important  feature  intimately  connected  with  some 
general  concepts  in  ODE  (Newman  &  Efroimsky  2003).  The  Hamilton-Jacobi  technique 
implicitly  fixes  the  gauge,  and  thus  leaves  the  internal  symmetry  heavily  veiled  (Efroimsky 
&  Goldreich  2003).  Below  we  shall  present  several  relevant  facts  and  formulae,  while  a  more 
comprehensive  treatment  of  the  matter  may  be  looked  up  in  (Efroimsky  &  Goldreich  2003, 
2004). 


As  well  known,  a  solution 

r  =  /  (Ci, ...,  C6,  t) 

(1) 

to  the  reduced  two-body  problem 

-  Gm  r 

r  +  —  -  =  ° 

(2) 

is  a  Keplerian  ellipse  or  hyperbola  parametrised  with  some  set  of  six  independent  orbital 
elements  C’i  which  are  constants  in  the  absence  of  disturbances. 

In  the  N-particle  setting  (or,  more  generally,  in  the  presence  of  whatever  disturbances) 
each  body  will  be  subject  to  some  perturbing  force  AF (r,  r,  t )  acting  at  it  from  bodies 
other  than  the  primary: 


r  + 


Gm  r 


AF 


(3) 


Solving  the  above  equation  of  motion  by  the  VOP  method  implies  the  use  of  (1)  as  an  ansatz, 


r  =  /(CW),  ..., 


(4) 


— f 

the  function  /  being  the  same  as  in  (1),  and  the  ’’constants”  Ci  now  being  endowed  with 
a  time  dependence  of  their  own.  Insertion  of  (4)  into  (3)  is  insufficient  for  determining  the 
six  functions  Ci(t )  .  because  the  vector  equation  (3)  comprises  only  three  scalar  equalities. 
To  furnish  a  solution  three  more  equations  are  needed.  Since  the  age  of  Lagrange  it  has  been 
advised  in  the  literature  to  employ  for  this  purpose  the  conditions  of  osculation, 


y-  Of  d( = 
i  dCi  dt 

imposition  whereof  ensures  that  the  physical  velocity, 

dr  _  Of  Of  dCj 

dt  dt  .  d  Ci  dt 

l  1 


(5) 

(6) 
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in  the  disturbed  case  is  expressed  by  the  same  function  g(Ci, 
two-body  configuration: 


g  = 


df 

dt 


... ,  C6,  t ) 


as  it  used  to  in  the 


(7) 


The  conditions  being  essentially  arbitrary,  their  choice  affects  only  the  mathematical  devel¬ 
opments  but  not  the  physical  orbit.  The  orbit’s  invariance  under  the  alternative  choices  of 
the  supplementary  constraints  reflects  the  gauge  freedom,  i.e.,  the  internal  symmetry  present 
in  this  problem.  In  the  language  of  pure  mathematics  this  is  a  fiber-bundle  structure.  In 
physicists’  terms  it  is  an  example  of  the  gauge  freedom.  Essentially,  this  is  merely  a  case  of 
ambiguous  parametrisation. 

The  Lagrange  constraint  and  the  law  of  motion,  (3),  with  ansatz  (4)  inserted  therein, 
yields  the  following  equation  of  the  elements’  evolution: 


E  [c„  c,\ 


dCj 

dt 


(8) 


[Cn  Cj\  being  the  matrix  of  Lagrange  brackets  introduced  as 


lr  r  1  df_  dg_  _  df_  dg_ 
[  n  ~  dCn  dCj  dCj  dCn 


(9) 


So  defined,  the  brackets  depend  neither  on  the  time  evolution  of  (7,  nor  on  the  choice  of 

— ¥  — ► 

supplementary  conditions,  but  solely  on  the  functional  form  of  /  (Cp...^  ,  t)  and  g  =  df  /dt. 
In  case  we  decide  to  relax  the  Lagrange  constraint  and  to  substitute  it  by 


v-  df  dCi  - 

X-acilt=  s(Cl . •'*> 


(10) 


$  being  some  arbitrary  function  of  time  and  parameters  Ci  (but,  for  the  reasons  of  sheer 
convenience,  not  of  time  derivatives  of  Ci),  then  instead  of  (8)  we  shall  get,  for  n  =  1,  ... ,  6  , 


E 


^  [Cn  Cj]  + 


df  d<f> 

dO,  d<7~ 


dCj 

dt 


dCn  dCn  dt  dCn 


(11) 


These  gauge-invariant  perturbation  equations  of  celestial  mechanics,  generated  by  direct 
application  of  the  VOP  method,  were  derived  in  (Efroimsky  2002b).  They  give  us  an  op¬ 
portunity  to  use,  in  an  arbitrary  gauge  <&,  the  Lagrange  brackets  (9)  defined  in  terms  of 

— ^ 

the  unperturbed  functions  /  and  g  .  Expressions  for  these  brackets  have  long  been  known, 
and  they  can  be  employed  to  write  down  the  planetary  equations  in  their  gauge-invariant 
form.  If  our  goal  were  to  arrive  to  the  customary  Lagrange  system  of  planetary  equations,  we 
would  now  stick  to  the  trivial  gauge  (5)  and  restrict  ourselves  to  conservative  forces  solely: 
AF  e*  dR(r ) / dr.  We  however  bear  in  mind  the  objective  of  exploring  the  interplay  of  two 
freedoms:  freedom  of  reference  frame  choice  and  that  of  gauge  fixing.  For  this  reason  not 
only  shall  we  leave  the  gauge  arbitrary  but  shall  also  reserve  for  our  disturbance  a  form 
general  enough  to  include  both  physical  and  inertial  forces,  the  latter  showing  themselves 
en  route  from  an  inertial  coordinate  system  to  an  accelerated  one.  An  example  at  hand 
is  a  satellite  orbiting  a  precessing  oblate  planet:  the  oblateness  will  give  birth  to  an  extra 
physical  force,  while  the  precession  will  result  in  inertial  forces  if  we  decide  to  describe  the 
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orbit  in  non-inertial  axes  associated  with  the  precessing  planet.  What  is  important  about 
the  inertial  inputs  is  that  they  are  not  only  position  but  also  velocity  dependent.  This  gives 
us  a  motivation  to  consider  velocity-dependent  disturbances.  Another  motivation  for  this 
comes  from  the  relativity  where  the  correction  to  Newton’  gravity  law  contains  velocities 
(Brumberg  1992) 

Expression  (11)  is  the  most  general  form  of  the  equations  for  the  orbital  elements,  in 
terms  of  the  disturbing  force.  To  get  the  generic  form  of  the  equations  in  terms  of  the 

Lagrangian  disturbance,  let  us  recall  that  the  (reduced)  two-body  problem  is  described  by 

^  ^  2  T  i, 

the  unit-mass  Lagrangian  CQ{ r,  ?,  t)  =  r  /2  —  U(r  ,  t),  canonical  momentum  p  =  r  , 

and  Hamiltonian  H0{ r,  p,  t)  =  p  2/2  +  U(r  ,  t)  ,  while  the  disturbed  setting  will  be 
described  by  the  perturbed  functions: 

•(  2 

£(r,  r,  t)  =  CQ  +  AC  =V—  -  U(r,  t)  +  A£(r,  r,  t)  ,  (12) 


and 


d  AC 

p  =  r  +  - —  , 

Or 


H( r,  p,  t)  —  pr  -  C  —  —  +  U  +  AH  , 


AH{r,  p,  t)  =  —  AC  — 


1  fdA  c\2 

2  i  dr  )  ’ 


(13) 


(14) 


AH  being  a  variation  of  the  functional  form:  AH  =  H( r,  p,  t)  —  7dc(r,  p,  t)  .  The 
Euler-Lagrange  equations  written  for  the  perturbed  Lagrangian  will  read: 


dU_ 

dr 


+  AF 


(15) 


the  new  term 


dAC  _  d_  (dAC\ 
8r  dt  \  dl  ) 


(16) 


being  the  disturbing  force.  We  see  that,  in  the  case  of  velocity-dependent  perturbations,  this 
force  is  equal  neither  to  the  gradient  of  the  Lagrangian’s  perturbation  nor  to  the  gradient 
of  negative  Hamiltonian’s  perturbation.  This  should  be  taken  into  account  when  comparing 
results  obtained  by  different  techniques.  For  example,  in  Goldreich  (1965)  the  word  “dis¬ 
turbing  function”  was  used  for  the  negative  perturbation  of  the  Hamiltonian.  The  gradient 
of  so  defined  disturbing  function  was  not  equal  to  the  disturbing  force. 

Substitution  of  (16)  into  (11)  then  yields  the  generic  form  of  the  equations  in  terms  of 
the  Lagrangian  disturbance.1  The  result, 


1  Direct  plugging  of  (16)  into  (11)  results  in 


^  ( ’"  ( dt 


df  dAL  df  d  W  dAL\ 
dCn  dr  8Cn  dt  \  +  / 


8g_ 

dCn 


$ 


5 


(17) 
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[Cn  Cj\  + 


of  o 

0Cn  OCj 


(  PAL 
V  dr 


dC3 

dt 


0 

OC~n 


A  L  +  - 
2 


/  rl  A/A2 

V  <9r  / 


dg  Of  0  dAL  3  \  0  AL\ 

0Cn  0Cn  dt  df  dCn  J  y  dr  ) 


not  only  reveals  the  convenience  of  the  generalised  Lagrange  gauge 


<£ 


d  AL 
dr 


(18) 


(which  reduces  to  =  0  in  the  case  of  velocity-independent  perturbations),  but  also  explic¬ 
itly  demonstrates  how  the  Hamiltonian  variation  comes  into  play:  it  is  easy  to  notice  that 
the  sum  in  square  brackets  is  equal  to  —  A7i(cont\  Substitution  of  the  general-type  force 
AF  into  (11)  or  of  the  Lagrangian  perturbation  AL  into  (17)  gives  us  the  gauge-invariant 
planetary  equations,  provided  we  know  the  expressions  for  elements  of  the  inverse  to  the 
Lagrange-bracket  matrix  [C(  Cj\  .  When  the  elements  Ci  are  chosen  as  the  Kepler  or 
Delaunay  variables,  (17)  entails  the  gauge-invariant  versions  of  the  Lagrange  and  Delaunay 
planetary  equations  (See  Appendix  1  to  (Efroimsky  &  Goldreich  2003).) 


2.2  Goldreich  (1965) 

In  1965  Goldreich  accomplished  a  ground-breaking  work  that  marked  the  beginning  of  studies 
of  the  Martian  satellite  dynamics.  He  started  out  with  two  major  assertions.  One  was  that 
the  Martian  satellites  had  either  been  formed  in  the  equatorial  plane  or  been  brought  therein 
very  long  ago.  The  second  was  that  Mars  has  experienced,  though  its  history,  a  uniform 
precession.  While  the  former  proved  to  be  almost  certainly  correct,2  the  validity  of  the  latter 
remains  model-dependent.  While  in  the  simplest  approximation  the  planetary  precession  is 

If  the  disturbance  were  to  depend  upon  positions  solely,  the  first  term  on  the  right-hand  side  of  the  above 
equation  would  be  equal  simply  to  8AL/dCn  .  In  general,  though,  we  should  rely  on  the  chain  rule 

d  AL  d  AL  df  d  AL  dr  d  AL  df  dAL  <9(g  +  4>) 
dCn  dr  dCn  +  df  dCn  dr  dCn  +  df  dCn 

insertion  whereof  into  the  preceding  formula  entails: 

Vfr  r]dCj  dAL  dAL  d$  (  d  fs  +  dAL  ) 

j  dt  dCn  gf  dCn  \ydCn  dt  dCn J  \  dr  ) 

To  arrive  herefrom  to  (17),  one  will  have  to  add  and  subtract  (l/2)8(8(AL)/8r)/8Cn  on  the  right-hand 
side.  Thereby  one  makes  the  gauge  function  $  appear  everywhere  in  the  company  of  d(AL)/8r  . 

2  Phobos  and  Deimos  give  every  appearance  of  being  captured  asteroids  of  the  carbonaceous  chondritic 
type  (Veverka  1977,  Pang  et  al.  1978,  Pollack  et  al.  1979,  Tolson  et  al.  1978).  Given  the  considerations 
summarised  in  detail  by  Murison  (1988),  it  seems  least  unlikely  that  the  Martian  satellites  were  formerly 
asteroids  that  were  captured  by  the  mechanism  of  gas  drag.  If  so,  this  must  have  occurred  early  in  the 
history  of  the  solar  system  while  the  gas  disk  was  substantial  enough  to  effect  a  capture.  At  that  stage  of 
planetary  formation,  the  spin  of  the  forming  planet  would  be  closely  aligned  with  the  orbital  plane  of  the 
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uniform,  a  more  involved  analysis,  carried  out  by  Ward  (1973,  1974),  Laskar  &  Robutel 
(1993),  and  Touma  &  Wisdom  (1994),  offers  evidence  of  strongly  non-uniform,  perhaps  even 
chaotic,  variations  of  the  Martian  obliquity  at  long  time  scales.  It  should  be  said,  though, 
that  the  analysis  presented  by  these  authors  was  mo  del- dependent.  In  particular,  it  was 
performed  in  the  approximation  of  the  planet  being  a  nondissipative  rigid  rotator.3  Since 
the  orbits  of  both  satellites  are  located  within  less  that  2°  from  the  equator,  the  two  said 
assertions  come  to  contradiction,  unless  there  exists  a  mechanism  constraining  satellite  orbits 
within  the  vicinity  of  the  primary’s  equator.  (Otherwise,  as  Goldreich  noted  in  his  paper, 
’’the  present  low  inclinations  of  these  satellites’  orbits  would  amount  to  an  unbelievable 
coincidence.”)  In  quest  of  such  a  mechanism,  Goldreich  (1965)  investigated  evolution  of  the 
Kepler  elements  of  a  satellite  in  a  reference  system  co-precessing  with  the  planet.  He  followed 
the  traditional  variation-of-parameters  (VOP)  scheme,  i.e.,  assumed  a  two-body  setting  as 
an  undisturbed  problem  and  then  treated  the  inertial  forces,  emerging  in  the  co-precessing 
frame,  as  perturbation  (along  with  the  other  perturbation  caused  by  the  oblateness  of  the 
planet).  Below  we  present  a  brief  summary  of  Goldreich’s  results,  with  only  minor  comments. 

Goldreich  began  with  applying  formulae  (12  -  16)  to  motion  in  a  coordinate  system 
attached  to  the  planet’s  centre  of  mass  and  precessing  (but  not  spinning)  with  the  planet. 
In  this  system,  the  equation  of  motion  includes  inertial  forces  and,  therefore,  reads: 

•4  dU  u  i.  .  _  ^  ,, 

or 

_  (19) 

planet’s  motion  about  the  Sun  i.e.,  the  obliquity  would  be  small  and  the  gas  disk  would  be  nearly  coplanar 
with  the  planetary  orbit.  Energetically,  a  capture  would  be  most  likely  to  be  equatorial.  This  is  most 
easily  seen  in  the  context  of  the  restricted  three- body  problem.  The  surfaces  of  zero  velocity  constrain  any 
reasonable  capture  to  occur  from  directions  near  the  inner  and  outer  collinear  Lagrange  points  (Szebehely 
1967,  Murison  1988),  which  lie  in  the  equatorial  plane.  Also,  a  somewhat  inclined  capture  would  quickly  be 
equatorialised  by  the  gas  disk.  If  the  capture  inclination  is  too  high,  the  orbital  energy  is  then  too  high  to 
allow  a  long  enough  temporary  capture,  and  the  object  would  hence  not  encounter  enough  drag  over  a  long 
enough  time  to  effect  a  permanent  capture  (Murison  1988).  Thus,  Phobos  and  Deimos  were  likely  (in  as 
much  as  we  can  even  use  that  term)  to  have  been  captured  into  near-equatorial  orbits. 

3  In  all  these  studies  the  planet  was  modelled  as  a  rigid  body.  Touma  &  Wisdom  (1994)  emphasised 
in  their  research  that  they  assumed  Mars  to  be  always  in  its  principal  rotation  state.  A  common  feature 
of  all  these  models  was  their  neglect  of  the  inelastic  relaxation  of  the  spin  axis’  precession.  Whether  the 
latter  simplification  is  physically  justified  is  yet  to  be  determined.  Inelastic  dissipation  is,  essentially,  the 
internal  friction  and  results  from  alternating  stresses  and  strains  emerging  in  a  wobbling  rotator  (Prendergast 
1958;  Burns  &  Safronov  1973;  Efroimsky  2001,  2002c;  Molina,  Moreno  &  Martinez-Lopez  2003;  Sharma, 
Burns  &  Hui  2004) .  It  is  well-known  that  the  lower  the  frequency  of  precession  the  slower  dissipation  of  the 
elastic  energy  associated  therewith.  This  circumstance  suggests  that  slow  obliquity  changes  will  result  in 
virtually  no  damping.  (The  same  circumstance  makes  one  think  that  the  fastest  modes  of  spin  precession 
are  well  dissipated.  It  was  for  this  reason  that  Touma  &  Wisdom  (1994)  accepted  the  model  of  Mars  always 
being  in  its  principal-rotation  state.)  The  real  physical  picture  is  far  more  complicated,  mainly  due  to  the 
essentially  nonlinear  nature  of  the  phenomenon.  Since  Mars  is  not  a  perfectly  symmetrical  oblate  spheroid 
but  is  a  triaxial  body,  its  free  rotation  in  a  spin  state,  which  is  even  slightly  different  from  the  principal 
one,  is  described  by  the  elliptic  functions  of  Jacobi  whose  decomposition  over  sines  and  cosines  contains 
an  infinite  multitude  of  harmonics.  From  here,  it  is  possible  to  demonstrate  that  in  the  course  of  rotation 
a  continuous  exchange  of  energy  among  these  harmonics  is  taking  place.  Due  to  this  phenomenon,  the 
alternating  stresses  associated  with  short-period  precession  are  always  getting  energy  from  stresses  caused 
by  long-period  changes  of  the  spin  axis.  This  way,  dissipation  of  the  short-time  motions  indirectly  leads 
to  damping  of  the  long-period  ones.  This  phenomenon,  called  “energy  cascade”  is  very  well  known  in  the 
theory  of  turbulence,  but  in  fact  it  is  generic  and  shows  itself  in  nonlinear  situations.  In  the  precessing-top 
context  it  has  not  been  explored  so  far,  and  we  do  not  know  how  effective  it  is  in  restraining  the  obliquity 
variations  caused  by  solar  torque. 
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p  x  (p  x  r)  , 


dUo 

dr 


d  AU 

dr 


2p  x  r  —  p  x  r  — 


clots  denoting  for  time  derivatives  in  the  co-precessing  frame  and  p  standing  for  the 
coordinate  system  angular  velocity  relative  to  an  inertial  frame.4  Here  the  physical  (i.e. , 
not  associated  with  inertial  forces)  potential  U (r)  consists  of  the  (reduced)  two-bocly  part 
U0(r)  =  —  GMr/r3  and  a  term  A U(r)  caused  by  the  planet’s  oblateness.  The  overall 
disturbing  force  on  the  right-hand  side  of  the  above  equation  is  generated,  according  to  (16), 
by 


AC  (r,  r,  t)  =  -  AU(r)  +  r-(pxr)  +  -  (p  x  ?)  •  (p  x  r) 

Zj 


Since  in  this  case 


then 


d  AC 

dr 


M  x  r  , 


^  dAC  u  ..  - 

p  =  r  +  - —  =  r  +  ^  x  r 

dr 

and,  therefore,  the  appropriate  Hamiltonian  variation  will  look: 


An  = 


„  1  fdAC\ 

“Hit 


—  —  AU  +  p  •  (p  x  r)  ]  =  AU  —  (r  x  p)  •  p  =  AU  —  J  ■  p  . 


(20) 


(21) 


(22) 


This  way,  A7d  becomes  expressed  through  quantities  defined  in  the  co-precessing  frame: 
the  satellite’s  orbital  momentum  vector  J  =  r  x  p  and  the  precession  rate  p  . 

Goldreich  employed  the  above  expression  in  the  role  of  a  disturbing  function  R  in  the 
planetary  equations: 

di  _  cos  i  9(  —  AH)  1  d  ( —  AH)  , 

dt  no?  (I  —  e2)1/2  sini  du  na2(  1  —  e2)1/2  sin*  dQ 


dn  =  1  d(-AH) 

dt  na2(  1  —  e2)1/2  sin  *  di 
where 

^  —  R  Roblate  I  Rinertial  (^6) 

4  Be  mindful  that  p,  ,  though  being  a  precession  rate  relative  to  an  inertial  frame,  is  a  vector  defined  in 
the  precessing  frame.  For  details  see  section  8.6  in  Marsden  and  Ratiu  (2003)  or  section  27  in  Arnold  (1989). 


8 


consists,  according  to  (22),  of  two  inputs:5 


R oblate  (  ^ )  —  A  f 


Gm  J2  p2 


1  —  3  sin2  i  sin2  (a;  +  u) 


(26) 


and 

Rinertial  =  i  ft  =  \J  G  m  CL  (1  -  e2)  W  fi  .  (27) 

Here  m  =  (mprirriary  +  msec(mdary)  •  The  mean  motion  is,  as  ever,  n  =  (G  m)1^2  a~3^2  , 
while  p  stands  for  the  mean  radius  of  the  primary,  v  denotes  the  true  anomaly,  and 


1  -  e2 

T  '  (1  - 

1  +  e  cos  u 


(28) 


is  the  instantaneous  orbital  radius.  In  the  right-hand  side  of  (27)  it  was  assumed  that  the 
angular  momentum  is  connected  with  the  orbital  elements  through  the  well-known  formula 

J  =  r  x  p  =  ^G  m  a  (1  —  e2)  w  (29) 


where 


w  =  xi  sin  i  sin  —  X2  sin  i  cos  O  +  X3  cos  i 


is  a  unit  vector  normal  to  the  instantaneous  ellipse,  expressed  through  unit  vectors  Xi,  x2,  x3 
associated  with  the  co-precessing  frame  xu  x%,  xs  (the  axes  x1  and  xs  lying  in  the  planet’s 
equatorial  plane).  The  afore  written  expression  for  w  evidently  yields: 


Rinertial  =  y  G  m  a  (l  —  c2)  (  pL\  sin  /  sin  H  —  p2  sin  i  cos  12  +  /i3  cos  i  ) 


(30) 


while  (26)  may  be,  in  the  first  approximation,  substituted  with  its  secular  part,  i.e.,  with  its 
average  over  the  orbital  period: 


(  Roblate  ) 


n 2  J2  2  3  cos  2i  —  1 

~T  P  (1  -  elf2  ' 


(31) 


the  averaging  having  been  carried  out  through  the  medium  of  formula  (109)  from  the  Ap¬ 
pendix,  with  (28)  inserted.  With  aid  of  (30)  and  (31),  the  planetary  equations  (24)  and  (25) 
will  simplify  to: 

d% 

—  =  —  Hi  cos  Cl  —  p2  sin  ,  (32) 


dt 


p cos i 
a)  (1  -  e2)2 


+  O 


Iffl 

■I 2  n 


(33) 


5  Our  formula  (26)  slightly  differs  from  the  one  employed  by  Goldreich  (1965),  because  here  we  use  the 
modern  definition  of  J2  : 


U 


E  T 

m= 2 


a  being  the  satellite’s  latitude  in  the  planet-related  coordinate  system.  The  coefficient  J  used  by  Goldreich 
(1965)  differs  from  our  J2  by  a  constant  factor:  J  =  (3/2)  J2  p2/r2  . 
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The  latter  results  in  the  well-known  node-precession  formula, 


n 


_  3  /  p\2  cos  i 

—  h n  2  It)  7i  2^  ^  ~  t°) 

i  \aj  (1  —  ez) 


(34) 


Its  insertion  into  the  former  entails 


i  =  -  —  cos  [  -  x  (t  -  t0)  +  D0]  +  —  sin  [  -  y  (t  -  ta)  +  QJ  +  ic 
X  v 


(35) 


where 


X 


cosi 

(1  -  e2)2 


(36) 


In  (Goldreich  1965),  equation  (35)  was  the  main  result,  its  derivation  being  valid  for  wobble 
which  is  slow  (  \p  |  J2  n  )  and  close  to  uniform  (\p,\/\p,\  ,/2  n ) . 

Despite  a  warning  issued  by  Goldreich  in  his  paper,  this  result  has  often  been  misinter¬ 
preted  and,  therefore,  misused  in  publications  devoted  to  satellites  and  rings  of  wobbling 
planets,  as  well  as  in  the  literature  on  orbits  about  tumbling  galaxies. 

In  (35)  i  stands  for  the  inclination  defined  in  co-precessing  axes  associated  with  the 
planet’s  equator,  and  therefore  (35)  clearly  demonstrates  that,  in  the  course  of  obliquity 
changes,  this  inclination  oscillates  about  zero,  with  no  secular  shift  accumulated.  Does  this 
necessarily  mean  that  the  satellite  orbit,  too,  oscillates  about  the  equatorial  plane,  without 
a  secular  deviation  therefrom?  Most  surprisingly,  the  answer  to  this  question  is  negative. 
The  reason  for  this  is  that  so  calculated  orbital  elements,  though  defined  in  the  co-precessing 
frame,  are  not  osculating  therein.  In  other  words,  in  the  frame  where  the  elements  are 
introduced,  the  instantaneous  ellipses  parametrised  by  these  elements  are  not  tangent  to  the 
physical  orbit  as  seen  in  this  frame. 

This  circumstance  was  emphasised  yet  by  Goldreich,  who  noticed  that  formula  (29) 
normally  (i.e. ,  when  employed  in  an  inertial  frame)  connects  the  osculating  elements  defined 
in  that  frame  with  the  angular  momentum  rxp  defined  in  the  same,  inertial,  frame  (i.e., 
with  rxr).  Since  in  the  above  calculation  the  frame  is  not  inertial  (and,  therefore,  the 
angular  momentum  is  different  from  rxr  but  is  equal  to  rxp  =  fx  (f  +  /Ixf)  ),  then 
the  orbital  elements  returned  by  (29)  cannot  be  osculating  in  this  frame.6  On  these  grounds 
Goldreich  warned  the  reader  of  the  peculiar  nature  of  the  elements  used  in  his  integration. 

To  this  we  would  add  that  it  is  not  at  all  evident  that  the  inertial-forces-caused  alteration 
of  the  planetary  equations  should  be  achieved  through  amending  the  disturbing  function  with 
the  momentum-dependent  variation  of  the  negative  Hamiltonian,  J  ■  p  .  While  the  common 
fallacy  identifies  the  disturbing  function  with  the  negative  Hamiltonian  perturbation,  in  real¬ 
ity  this  rule-of-thumb  works  (and  yields  elements  that  are  osculating)  only  for  disturbances 
dependent  solely  upon  positions,  not  upon  velocities  (i.e.,  for  Hamiltonian  perturbations 


6  Were  these  elements  osculating  in  the  frame  wherein  they  had  been  defined,  then  formula  (29)  would 
read:  r  x  r  =  \/Gm  a  (1  —  e2)  w  ,  i.e.,  would  connect  the  elements  with  the  velocity  in  that  frame. 
In  reality,  though,  it  reads:  rxp  =  y/Gm  a  (1  —  e2)  w  ,  i.e.,  connects  the  elements  with  the 
momentum  p  =  r  +  p,  x  r  which  happens  to  coincide  with  the  satellite’s  velocity  relative  to  the  inertial 
axes.  This  situation  was  formulated  by  Goldreich  in  the  following  terms:  the  orbital  elements  emerging 
in  the  above  derivation  are  defined  in  the  co-precessing  frame  but  are  osculating  in  the  inertial  one.  This 
illustrative  metaphor  should  not,  though,  be  overplayed:  the  fact  that  the  elements  emerging  in  Goldreich’s 
computation  return  the  inertial- frame-related  velocity  does  not  mean  that  this  inclination  may  be  interpreted 
as  that  relative  to  the  invariable  plane.  (The  elements  were  introduced  in  the  co-precessing  frame!) 
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dependent  only  upon  coordinates,  but  not  upon  momenta).  Ours  is  not  that  case  and, 
therefore,  more  alterations  in  the  planetary  equations  are  needed  to  account  for  the  frame 
precession,  if  we  wish  these  equations  to  render  osculating  elements.  However,  if  one  neglects 
this  circumstance  and  simply  amends  the  disturbing  function  with  J  •  jl ,  then  the  planetary 
equations  will  give  some  elements  different  from  the  osculating  ones.  It  will  then  become  an 
interesting  question  as  to  whether  such  elements  will  or  will  not  coincide  with  those  rendered 
by  (28)  when  this  formula  is  used  in  non-inertial  frames. 

All  these  subtle  issues  get  untangled  in  the  framework  of  the  gauge  formalism.  Appli¬ 
cation  of  this  formalism  to  motions  in  non-inertial  frames  of  references  was  presented  in 
(Efroimsky  &  Goldreich  2004).  The  main  results  proven  there  are  the  following. 

1.  If  one  attempts  to  account  for  the  inertial  forces  by  simply  adding  the  term  J  •  fl 
to  the  disturbing  function,  with  no  other  alterations  made  in  the  planetary  equations,  then 
these  equations  indeed  do  render  quantities  that  may  be  interpreted  as  some  orbital  elements 
(i.e.,  as  parameters  of  some  instantaneous  conics).  These  elements  are  NOT  osculating  and, 
therefore,  the  instantaneous  conics  parametrised  by  there  elements  are  not  tangent  to  the 
physical  orbit.  Hence,  these  elements  cannot,  generally,  be  attributed  a  direct  physical  in¬ 
terpretation,7  except  in  the  situations  when  their  deviation  from  the  osculating  elements 
remains  sufficiently  small. 

2.  By  a  remarkable  coincidence,  these  non-osculating  elements  turned  out  to  be  identical 
with  those  emerging  in  formula  (29).  This  coincidence  was  implicitly  taken  for  granted  by 
Goldreich  (1965),  which  reveals  his  truly  incredible  scientific  intuition. 

3.  To  build  up  a  system  of  planetary  equations  that  render  osculating  elements  of  the 
orbit  as  seen  in  the  co-precessing  coordinate  system,  one  has  not  only  to  add  J  •  fl  to  the 
disturbing  function,  but  also  to  amend  each  of  these  equations  with  several  extra  terms. 
Some  of  those  terms  are  of  order  ( \fl  |/( J2  n)  )2  ,  some  others  are  of  order  \fl\/(\fl\  J2  n) . 
Most  importantly,  some  terms  are  of  the  first  order  in  the  precession-caused  perturbation 
| fl  |/(J2  n)  ,  which  means  right  away  that  the  non-osculating  elements  used  by  Goldreich 
(1965)  differ  from  the  osculating  ones  already  in  the  first  order.  While  more  comprehensive 
account  on  this  topic,  with  the  resulting  equations,  will  be  offered  in  the  end  of  this  section, 
here  we  shall  touch  upon  only  one  question  which  gets  immediately  raised  by  the  presence 
of  such  first-order  differences.  This  question  is:  what  are  the  averages  of  these  differences? 
Stated  alternatively,  do  the  secular  components  of  the  said  non-osculating  elements  differ 
considerably  from  those  of  the  osculating  ones?  Goldreich  (1965)  stated,  without  a  proof, 
that  the  secular  components  differ  only  in  high  orders  over  the  velocity- dependent  part  of 
the  perturbation.  In  our  paper  we  shall  probe  the  limits  for  this  assumption. 


2.3  Brumberg  and  Kinoshita 

A  development,  part  of  which  was  similar  to  that  of  Goldreich  (1965),  was  independently 
carried  out  by  Brumberg,  Evdokimova  &  Kochina  (1971)  who  studied  orbits  of  artificial  lunar 

7When  the  orbit  evolution  is  sufficiently  slow,  the  observer  can  attribute  some  physical  meaning  to 
elements  of  the  osculating  conic.  For  example,  whenever  an  observer  talks  about  the  inclination  or  the 
eccentricity  of  a  perturbed  orbit,  he  naturally  implies  those  of  the  osculating  ellipse  or  hyperbola. 
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satellites  in  a  coordinate  system  co-precessing  with  the  Moon.  In  that  article,  too,  the  non¬ 
osculating  nature  of  the  resulting  orbital  variables  did  not  go  unnoticed.  The  authors  called 
these  variables  “contact  elements”  and  stated  (though  never  proved)  that  these  variables 
return  not  the  correct  value  of  the  velocity  but  that  of  the  momentum.  Later,  one  of  these 
authors  rightly  noted  in  his  book  (Brumberg  1992)  that  the  contact  elements  differ  from  the 
osculating  ones  already  in  the  first  order  over  the  velocity-dependent  part  of  the  perturbation. 
In  subsection  1.1.3  of  that  book,  he  unsuccessfully  tried  to  derive  analytical  transformations 
interconnecting  these  sets  of  variables.8 

A  similar  attempt  was  undertaken  in  a  very  interesting  article  by  Ashby  &  Allison  (1993). 
Though  the  authors  succeeded  in  many  other  points,  their  attempt  to  derive  formulae  for 
such  a  gauge  transformation  was  not  successful.9 

The  setting,  considered  by  Goldreich  (1965)  in  the  context  of  Martian  satellites  and  by 
Brumberg  et  al  (1971)  in  the  context  of  circumlunar  orbits,  later  emerged  in  the  article  by 
Kinoshita  (1993)  who  addressed  the  satellites  of  Uranus. 

Kinoshita’s  treatment  of  the  problem  was  based  on  the  following  mathematical  construc¬ 
tion.  Denote  satellite’s  positions  and  velocities  in  the  inertial  and  in  the  co-precessing  axes 
with  {r  ' ,  v  '  }  and  with  {  r  ,  v}  ,  correspondingly.10  Interconnection  between  them  will 
be  implemented  by  an  orthogonal  matrix  A  , 

r  =  if'  ,  v  =  r  =  Ar  '  +  Ar  '  =  A  A~l  r  I  A  v  '  ■-  -  —  p  x  r  +  A  p  '  ,  (37) 

p  being  the  precession  rate  as  seen  in  the  co-precessing  coordinate  system,  and  the  inertial 
velocity  v  '  being  identical  to  the  inertial  momentum  p  '  .  Kinoshita  suggested  to 
interpret  this  interconnection  as  a  canonical  transformation  between  variables  {  r  ' ,  p  '  } 
and  {  r  ,  p  }  ,  implemented  by  generating  function 

F2  =  p  ■  Ar'  =  (yAA p^j  •  r'  .  (38) 

This  choice  of  generating  function  rightly  yields 

w =  Ar'  ’  (39) 

8  Contrary  to  the  author’s  statement,  formula  (1.1.41)  in  (Brumberg  1992)  is  not  rigorous,  but  is  valid 
only  in  the  first  order.  (To  make  it  rigorous,  one  should  substitute  everywhere,  except  in  the  denominators, 
r  with  r  —  dR/dr  .)  Besides,  the  author  did  not  demonstrate  his  derivation  of  formula  (1.1.43),  for  Ma  , 
from  (1.1.42).  (In  Brumberg’s  book  the  mean  anomaly  is  denoted  with  /,  not  with  M  .  ) 

Most  importantly,  the  qualitative  reasoning  presented  by  the  author  in  the  paragraph  preceding  formula 
(1.1.43)  is  un-rigorous  and  essentially  incorrect.  The  cause  of  this  is  that  the  author  compares  the  planetary 
equations  for  contact  elements,  written  in  a  precessing  frame,  with  the  equations  for  osculating  ones,  written 
in  an  inertial  frame,  instead  of  comparing  two  such  systems  (for  contact  and  for  osculating  elements)  both 
of  which  are  written  in  a  precessing  frame.  This  makes  a  big  difference  because,  as  we  already  explained 
above,  transition  to  a  precessing  frame  does  not  simply  mean  addition  of  an  extra  term  to  the  disturbing 
function. 

Despite  all  these  mathematical  irregularities,  the  averaged  system  of  planetary  equations  (1.1.44),  “de¬ 
rived”  by  Brumberg  for  the  first-order  secular  perturbations,  turns  out  to  be  correct  in  the  limit  of  uniform 
precession.  Just  as  in  the  preceding  subsection  we  had  a  reason  to  praise  the  unusual  intuition  of  Goldreich, 
so  here  we  have  to  pay  tribute  to  the  exellent  intuition  of  Brumberg,  intuition  which  superseded  his  flawed 
mathematics. 

9  To  carry  out  the  gauge  transformation,  the  authors  used  a  set  of  intermediate  variables  {Q(0) ,  Pk  (o)} 
which  were  canonical  and,  at  the  same  time,  osculating.  As  follows  from  the  theorem  proven  by  Efroimsky 
&  Goldreich  (2003),  these  variables  are  nonexistent  when  the  perturbation  depends  upon  velocities. 

10  We  use  notations  opposite  to  those  in  Kinoshita  (1993),  in  order  to  conform  with  the  notations  of 
Goldreich  (1965). 
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(40) 


while  the  interconnection  between  momenta  will  look: 

P'  =  |^7  =  AT p  ,  i.e.,  p  =  (iT)  1  p  '  =  Ap'  , 
whence  p  =  v  +  P  x  r  .  The  Hamiltonian  in  precessing  axes  will  read 

H(r ,  p)  =  H{inert)  (r  ' ,  p')  +  =  //(merf)  (r ',  p')  +  p  Ar' 

m  H{inert)  (r  ',  p  ')  -  p  ■  (p  x  r)  =  H{inert)  (r  ',  p  ')  -  (r  x  p)  ■  p  .  (41) 

The  Hamiltonian  perturbation,  caused  by  the  inertial  forces,  is  —  (r  x  p)  p  —  J  x  /I , 
vector  J  being  the  orbital  momentum  as  seen  in  the  co-precessing  frame.  Comparing  this 
with  (22),  we  see  that  employment  of  the  above  canonical  transformation  is  but  another 
method  of  stepping  on  the  same  rake.  In  distinction  from  Goldreich  (1965)  and  Brumberg 
et  al  (1971),  Kinoshita  in  his  paper  did  not  notice  that  he  was  working  with  non-osculating 
elements. 

The  problem  with  Kinoshita’s  treatment  is  that  the  condition  of  canonicity  in  some  sit¬ 
uations  comes  into  contradiction  with  the  osculation  condition.  In  other  words,  canonicity 
sometimes  implicitly  contains  a  constraint  which  sometimes  is  different  from  the  Lagrange 
constraint  (5).  This  issue  was  comprehensively  elucidated  in  the  work  Efroimsky  &  Goldreich 
(2003).  The  authors  began  with  the  reduced  two-body  setting  and  thoroughly  re-examined 
the  Hamilton- Jacobi  procedure,  which  leads  one  from  the  spherical  coordinates  and  the  cor¬ 
responding  canonical  momenta  to  the  set  of  Delaunay  variables.  While  in  the  undisturbed 
two-body  case  this  procedure  yields  the  Delaunay  variables  which  are  trivially  osculating 
(and  parametrise  a  fixed  Keplerian  ellipse  or  hyperbola),  in  the  perturbed  case  the  situa¬ 
tion  becomes  more  involved.  According  to  the  Theorem  proven  in  that  paper,  the  resulting 
Delaunay  elements  are  osculating  (and  parametrise  a  conic  tangent  to  the  perturbed  tra¬ 
jectory)  if  the  Hamiltonian  perturbation  depends  solely  upon  positions,  not  upon  momenta 
(or,  the  same:  if  the  Lagrangian  perturbation  depends  upon  positions  but  not  upon  veloci¬ 
ties).  Otherwise,  the  Delaunay  elements  turn  out  to  be  non-osculating  (and  parametrise  the 
physical  trajectory  with  a  sequence  of  non-tangent  conics).  As  one  can  see  from  the  above 
equation,  the  Hamiltonian  perturbation,  caused  by  the  inertial  forces,  depends  upon  the 
momentum,  and  this  circumstance  indicates  the  problem.  This  trap,  in  which  many  have 
fallen,  is  of  a  special  importance  in  the  General  Relativity,  because  the  relativist  corrections 
to  the  equations  of  motion  are  velocity-dependent.11 


3  Planetary  equations 

In  this  section  we  shall  briefly  spell  some  results  obtained  in  Efroimsky  &;  Goldreich  (2004) 
and  shall  use  these  results  to  derive  the  Lagrange-type  planetary  equations  (59  -  64)  for 

11  In  an  interesting  article  (Chernoivan  &  Mamaev  1999),  the  authors  addressed  the  two-body  problem  on 
a  curved  background.  The  curvature  entailed  a  velocity-dependent  relativist  correction,  which  was  treated 
as  a  perturbation.  After  carrying  out  the  Hamilton- Jacobi  development,  the  authors  arrived  to  canonical 
variables  analogous  to  the  Delaunay  elements.  Orbit  integration  in  terms  of  these  variables  would  be  as 
correct  as  in  terms  of  any  others.  The  problems  began  when  the  authors  used  these  elements  to  arrive  to 
some  conclusions  regarding  the  perihelion  precession.  Those  conclusions  need  to  be  reconsidered,  because 
they  were  rendered  on  the  basis  of  Delaunay  elements  that  were  non-osculating.  Similar  comments  may  be 
made  about  the  work  by  Richardson  &  Kelly  (1988)  who  addressed,  using  the  Hamiltonian  language,  the 
two-body  problem  in  the  post-Newtonian  approximation. 
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osculating  elements  in  a  coordinate  system  co-precessing  with  an  oblate  primary. 


3.1  Planetary  equations  for  contact  elements 

Above,  in  subsection  2.1,  we  provided  a  very  short  account  of  the  gauge  formalism.  Ex¬ 
pression  (17),  presented  there,  is  the  most  general  form  of  the  planetary  equations  for  an 
arbitrary  set  of  six  independent  orbital  elements,  written  in  terms  of  arbitrary  disturbance 
of  the  Lagrangian. 

When  the  elements  (7*  are  chosen  to  be  the  Keplerian  or  Delaunay  sets  of  variables, 
we  arrive  to  the  gauge- invariant  versions  of  the  Lagrange  or  Delaunay  planetary  equations, 
correspondingly.  They  are  written  down  in  Appendix  to  the  paper  (Efroimsky  &  Goldreich 
2003).  Interplay  between  the  gauge  freedom  and  the  freedom  of  frame  choice  is  explained 
at  length  in  Section  3  of  the  article  (Efroimsky  &  Goldreich  2004)  which  addressed  orbits 
about  a  precessing  planet.  It  is  demonstrated  in  that  work  that,  if  one  chooses  to  describe 
the  motion  in  terms  of  the  non-osculating  elements  that  were  introduced  in  a  co-precessing 
frame  and  were  defined  in  the  generalised  Lagrange  gauge12  (18),  then  the  corresponding 
Hamiltonian  perturbation  will  read: 


^j^icont) 


Roblate(f)  +M'(/Xg) 


(42) 


while  the  planetary  equations  (17)  acquire  the  form 

dCi  _  d  (  -  A H^nt)  ) 


or: 


[Cr  Q 


da 


dt 


d 


'C'C‘)lF  =  9C 


9Cr 


Roblateif)  +  M  •  (/  X  g) 


(43) 


where  /  and  g  stand  for  the  un-disturbed  (two-body)  functional  expressions  of  the  position 
and  velocity  via  the  time  and  the  chosen  set  of  orbital  elements: 


/  ((A, ...,  Ce,  t ) 


(44) 


d 


v  =  g(C'1,...,C'6,  t)  =  —f{Cu...,C6,t) 


— f 

(so  that  /  (Ci(f), ...,  Ce(t ),  t)  and  g  (C'i(t), ...,  Ce(t),  t )  become  the  ansatz  for  solving  the 
disturbed  problem).  For  Relate  in  (42)  -  (43),  one  can  employ,  dependent  upon  the  desired 
degree  of  rigour,  either  the  exact  expression  (26)  or  its  orbital  average  (31). 

In  the  generalised  Lagrange  gauge  (18)  the  canonical  momentum  becomes: 


^  dA  C  ,  -j  dAC 

p  =  r  +  - -  g  <I>  - — 

dr  dr 

which  means  that  its  functional  dependence  upon  the  time  and  the  chosen  set  of  orbital 
elements  is  the  same  as  it  used  to  be  in  the  unperturbed  case  where  both  the  velocity  and 

12  It  is  an  absolutely  crucial  point  that  choice  of  a  gauge  and  choice  of  a  reference  frame  are  two  totally 
independent  procedures.  In  each  frame  one  has  an  opportunity  to  choose  among  an  infinite  variety  of  gauges. 
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the  momentum  were  simply  equal  to  g  (Cb, ...,  C6,  I)  .  This  also  means  that 
coincides  with  Goldreich’s  AH  given  by  (26). 

We  see  that  the  generalised  Lagrange  gauge  (18)  singles  out  the  same  set  of  non-osculating 
elements  which  showed  up  in  the  studies  by  Goldreich  (1965)  and  Brumberg,  Evdokimova 
&  Kochina  (1971),  -  the  set  of  “contact  elements.”  This  is  why  in  (42)  the  Hamiltonian 
perturbation  was  written  with  the  superscript  “cont.” 

In  (43)  the  Lagrange-bracket  matrix  is  defined  in  the  unperturbed  two-body  fashion 
(9)  and  can,  therefore,  be  trivially  inverted.  Hence,  when  the  elements  are  chosen  as  the 
Keplerian  ones,  the  appropriate  equations  look  like  the  customary  Lagrange-type  equations 
(i.e.,  like  (24)  and  (25)  above),  with  the  disturbing  function  given  by  (26)  or,  the  same,  by 
(42): 


da  _  _2_ 

dt  na  dM0 


(46) 


de  1-e2  '>  (  -  AH1""'"  )  (i  _  e2)i/2  8  (  - 
dt  na2  e  dM0  n a2  e  du 


(47) 


dw  _  -cos i  d(-AH^)  (i_e2)i/2  d(-AH(^) 

dt  na2(  1  —  e2)1/2  sin  i  di  na2e  de 


(48) 


di  _  cos i  _ 1  d(-AH^) 

dt  na2(  1  —  e2)1/2  sini  du  na2(  1  —  e2)1/2  sini  dtt 


dtt  _  1  d(-AH^) 

dt  na2(  1  —  e2)1/2  sini  di 


dM0 

dt 


1  _  e2  d  (  -  A H{cont))  2  <9  ( -  A H{cont)) 
na2  e  de  na  da 


(51) 
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3.2  Planetary  equations  for  osculating  elements 

When  one  introduces  elements  in  the  precessing  frame  and  also  demands  that  they  osculate 

— ^ 

in  this  frame  (i.e.,  makes  them  obey  the  Lagrange  gauge  $  =  0)  then  the  Hamiltonian 
variation  reads:13 

A  H^  =  -  [  Roblate  iy)  -f  £  •  (/  X  g)  +  (M  x  f)  .(/IX/)]  ,  (52) 

while  equation  ( 17 )  becomes: 


To  ease  the  comparison  of  this  equation  with  (43),  it  is  convenient  to  split  the  expression  for 
in  (52)  into  two  parts: 

An(cont)  =  _  J  t)  +  ft-  (/ X  g)  ]  (54) 

and 

-  (M  x  /)  •  (jZ  x  /)  ,  (55) 

and  then  to  group  the  latter  part  with  the  last  term  on  the  right-hand  side  of  (53): 


One  other  option  is  to  fully  absorb  the  0(\f2\2)  term  into  AH  ,  i.e.,  to  introduce  the 
effective  “Hamiltonian” 

AH(eff)  =  -  Robiate{v)  +/*'(/xg)  +  ^  {ft  x  /)  •  (p,  x  /)  (57) 

and  to  write  down  the  equations  like  this: 


\cn  ct\  ^ 


d  AH(eff) 

dCn 


x  g  ~  f  x 


•  (58) 


13  Just  as  A Ji^cont)  in  (42),  this  Hamiltonian  variation  is  still  equal  to  —  t)  +  ft  •  Jj  = 

-  R(f,  t)  +  p  (f  x  p)  .  However,  the  canonical  momentum  now  is  different  from  g  and  reads  as: 
P  =  g  +  (/*x  f)  ■ 
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For  Ci  being  chosen  as  the  Keplerian  elements,  inversion  of  the  Lagrange  brackets  will 
yield  the  following  Lagrange-type  system: 


da  2 

dt  na 


S(-A«Wfl)  /_x  8f 


dM0 


d  M0 


(59) 


de  1  —  e2 


dt  n  a 2  e 


d  M0 


dM0 


(60) 


(f  -  e2)1/2 

na2  e 


d(-A  ^  (df 


du 


M'hrxs-/X^r 

<9a;  <9u; 


i  v  df 


da; 

dt 


—  cos  i 


na2(  1  —  e2)1/2  sin  i 


'a(-Awwfl)  ^  (df  ^  -  9g 

- a -  +  "■  a  xe-f*di 


-  M  •  /  X 


5/ 

<9i 


(61) 


(1  -  e2)1/2 
na2  e 


<9  ( -  AH{eff)^ 


de 


_  .  df  _  ?  d 

+  (‘lSxg-/x-8e 


.  <9/ 

"  l/xs 


di  cos  i 

dt  na2  (1  —  e2)1/2  sin  i 


a(-AWw/>)  (df  ^  gg\  (  a/ 

— a: —  +  '*•  Ujxg-/xa;  -'*•  /xs: 


(62) 


na2  (1  —  e2)1/2  sin  i 


8(-A «<•">)  _  /a/  ,  -  8g\  ^  (T  df 

— m —  +  "  Mxg^/xM  -'*•  /XS2 


dCl 

dt 


na 2  (1  —  e2)1/ 2  sin  i 


d  (  -  A H(eff)  ) 


di 


__  df 

+  f,|w xg 


/X  f 


,  v  <9/ 

"■|/x  a 


,  (63) 
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(64) 


dMa 

dt 


1  -  e2 
no?  e 


d  (  -  A n{eff)) 

de 


+  P  ' 


2 

'  d  (  -  AH(eff)) 

n  a 

da 

P 


terms  /2  •  (  ( df  jdM0 )  x  g  —  ( dg/dM0 )  x  /  )  being  omitted  in  (59  -  60),  because  these 
terms  vanish  identically  (see  the  Appendix).  In  equations  (58  -  64),  AH is  given  by 
(57).  In  a  rigorous  analysis,  with  the  true-anomaly  dependence  of  all  terms  in  (59  -  64)  taken 
into  account,  the  Hamiltonian  will  look,  according  to  (26  -  27)  and  (30): 


AH(eff) 


(iA 


1  / 

/  -*■  iA 

2  (mx/) 

•p/j 

Gm  J2  p2 


1  —  3  sin2  i  sin2(tu  +  v)  +  \J  G  m  a  (1  —  e2)  w  ■  p  +  -  [p  x  f^j  ■  (ji  x 


G  m  J2  p 2 


1  +  e  cos v 
1  -  e2 


1  —  3  sin2  i  sin2  (a;  +  u) 


+ 


\jGm  a  (1  —  e2) 


(  pL\  sin  i  sin  12 


p2  sin  i  cos  12  +  cos  i  )  + 


(65) 


+  \  (m  x  f)  ■  (p  x  j)  . 

Otherwise,  when  all  terms  in  the  right-hand  sides  of  (59  -  64)  are  substituted  with  their 
orbital  averages  (denoted  with  the  (...)  symbol),  the  Hamiltonian  will  be,  due  to  (31): 


AH{eff) 


{ Robiate )  +  P  ■  (/x  g)  +  ^  (  (p  x  f)  ■  (p  x  /)  ) 


(66) 


G  m  J2  p2  3  cos2 1  —  1  rz - 77 - —  .  .  .  _  .  .  _ 

-  —7  - +  \  Gma  1  —  ez)  p,i  sinj  sml2  —  u2  snu  cosiz  +  u3  cos?)  + 

4  a3  (1  _  e2)3/2  v  v  7  v 


+  \  (  (m  x  f)  ■  (P  x  f)  ) 
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The  expression  for  /  x  g  is  true-anomaly-independent  and,  therefore,  does  not  need  to 
be  bracketed  with  the  averaging  symbols  (...}.  The  expression  for  (^p  x  ■  {p  x  f\ 
through  the  orbital  elements  is  too  cumbersome,  and  here  we  do  not  write  it  down  explicitly. 
(See  formula  (186)  in  the  Appendix.)  When  we  permit  ourselves  to  neglect  the  0( \P\2) 
inputs,  all  three  Hamiltonians,  A7Posc^  ,  A7i(conl'>  ,  and  A7PefR  coincide: 


a n(eff)  »  a n(osc)  a n(cont) 


(Roblate)  +  P  (/  X  g) 


(67) 


G  m  J2  p2  3  cos2  i  —  1 

4  a3  (i  _  e2)3/2 

Two  important  issues  should  be  dwelt  upon  at  this  point. 

First,  we  would  remind  that  the  function  AH(cont 1  ,  given  by  expression  (42),  yields 
the  correct  functional  form  of  the  Hamiltonian  only  in  the  case  we  express  the  Hamiltonian 
through  the  contact  elements  (and  calculate  these  through  (43)  or  (46  -  51)).  In  the  currently 
considered  case  of  osculating  elements,  this  A7d -con^  is  NOT  the  correct  expression  for  the 
Hamiltonian.  The  correct  functional  dependence  of  the  Hamiltonian  upon  the  osculating 
elements,  A 7 Posc">  ,  is  given  by  formula  (52).  Though  in  this  dynamical  problem  the  Hamil¬ 
tonian  is  unique,  its  functional  dependencies  through  the  contact  and  through  the  osculating 
elements  differ  from  one  another,  one  being  A7d'co”0  as  in  (42),  another  being  AH^0^  as 
in  (52).  As  for  the  function  AH^^  rendered  by  (57),  it  is  not  really  a  Hamiltonian,  but 
is  simply  a  convenient  mathematical  entity.  In  the  approximation,  where  0(\p\2)  terms 
are  neglected,  there  is  no  difference  between  these  three  functions.  Despite  this,  the  0(\p\) 
and  0(\p\)  terms  do  stay  in  equations  (58  -  64)  for  osculating  elements. 

Second,  we  would  comment  on  our  use  of  expressions  (29  -  30)  in  our  derivation  of  (65 
-  67).  Above,  in  subsections  2.2  and  3.1,  the  use  of  (27)  and  (30)  was  based  on  formula 
(29)  which  interconnected  J  r  x  p  =  f  x  (r  +  p  x  with  contact  elements 
a  ,  e  ,  i  ,  and  Q  .  As  demonstrated  in  Efroimsky  &  Goldreich  (2004),  in  that  frame 
r  =  g  —  p  x  /  .  Hence,  formula  (29)  interconnects  the  contact  elements  with 

— >  — f  a 

J  =  /  x  P  ■  In  the  present  subsection,  we  use  formula  (29)  in  its  usual  capacity,  i.e.,  to 
interconnect  J  r  x  r  f  x  g  with  osculating  elements  a  ,  e  ,  i  ,  and  Q  . 
It  may  seem  confusing  that,  though  in  both  cases  this  formula  can  be  written  down  in  the 
same  way, 


+  \JCma  (1  —  e2)  ( //,]  sin i  sin  D  —  p2  sin i  cos +  p3  cos i ) 


f  x  g  =  \JCm  a  (1  —  e2)  w  ,  (68) 

its  meaning  is  so  different.  The  clue  to  understanding  this  difference  lies  in  the  fact  that  in 
one  case  r  g  —  p  x  f  (and,  therefore,  the  elements  are  contact),  while  in  the  other 
case  r  =  g  (which  makes  the  elements  osculating).  For  more  details  see  Efroimsky  & 
Goldreich  (2004). 
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4  Comparison  of  the  orbital  calculations,  performed  in 
terms  of  the  contact  elements,  with  those  performed 
in  terms  of  the  osculating  elements.  The  simplest 
approximation. 

As  explained  in  Section  2,  it  follows  from  equations  (23  -  24)  that  the  initially  small  incli¬ 
nation  remains  so  in  the  course  of  the  oblate  primary’s  precession.  Whether  this  famous 
result  may  be  interpreted  as  keeping  of  satellites  in  the  near-equatorial  zone  of  a  precessing 
planet  will  depend  upon  how  well  the  non-osculating  (contact)  inclination  emerging  in  (23  - 
36)  approximates  the  physical,  osculating,  inclination  rendered  by  (59  -  64). 


4.1  The  averaged  planetary  equations. 

Comparing  equations  (59  -  64)  for  osculating  elements  with  equations  (46  -  51)  for  contact 
elements,  we  immediately  see  that  they  differ  already  in  the  first  order  over  the  precession 
rate  //  and,  therefore,  the  values  of  the  contact  elements  will  differ  from  those  of  their 
osculating  counterparts  in  the  first  order,  too.  A  thorough  investigation  of  this  difference 
would  demand  numerical  implementation  of  both  systems  and  would  be  extremely  time- 
consuming.  Meanwhile,  we  can  get  some  preliminary  estimates  by  asking  the  following, 
simplified,  question:  how  will  the  secular,  i.e.,  averaged  over  an  orbital  period,  components 
of  the  contact  and  orbital  elements  differ  from  one  another?  To  answer  this  question,  we 
perform  the  following  approximations: 

H  In  equations  (59  -  64),  we  substitute  both  the  Robiate  term  in  the  Hamiltonian  and 
the  //-dependent  terms  with  their  averages  (so  that,  for  example,  the  Robiate  term  will  be 
now  substituted  with  ( Robiate )  expressed  via  (31)  ). 

2.  We  neglect  the  terms  of  order  ft 2 .  This  way,  we  restrict  the  length  of  time  scales 
involved.  (Over  sufficiently  long  times  even  small  terms  may  accumulate  to  a  noticeable 
secular  correction.)  However,  we  can  now  benefit  from  the  approximate  equality  (67). 

So  truncated  and  averaged  system  of  Lagrange-type  equations  will  read: 


da  2 
dt  na 


(69) 


de  1  —  e2 

dt  n  a 2  e 


) 


(1  -  e2)1/2 
na2  e 


(  R 


(70) 
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duj  —  cos  % 

dt  na2(l  —  e2)1/2  sin  i 


a(-A «<«">)  (af  _  dg\  j  9f 

- ai -  +  <MaIxg  -fxel  >  ~  (f  x  si |} 


(71) 


+ 


(1  -  e2)1/2 
na 2  e 


di  _  cos  i 
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2 
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the  Hamiltonian  here  being  approximated  by  (67),  and  the  angular  brackets  signifying  or¬ 
bital  averaging.  In  equations  (69),  (70)  and  (72)  we  took  into  account  that  the  averaged 
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Hamiltonian  (67)  bears  no  dependence  upon  M0  and  u>  .  (This,  though,  will  not  be  the 
case  for  the  exact,  ^-dependent,,  Hamiltonian  given  by  (52)  and  (26)  !  ) 

Calculation  of  (  ( df/dCj )  x  g  —  (dg/dCj)  x  /  j  and  —  p,  (f  x  dg/dCj  )  takes 
pages  of  algebra.  A  short  synopsis  of  this  calculation  is  offered  in  the  Appendix.  Here  follows 
the  outcome: 
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Hi ,  H2  ,  A* 3  being  the  Cartesian  components  of  the  precession  rate,  as  seen  in  the  co- 
precessing  frame,  and  [j,±  being  the  component  of  the  precession  rate,  aimed  in  the  direction 
of  the  angular  momentum;  it  is  given  by  (117). 

It  may  seem  strange  that  the  right-hand  side  of  (76)  does  not  vanish  in  the  limit  of 
e  — »  0  .  The  absurdity  of  this  will  be  easily  redeemed  by  the  fact  that  this  term  shows  up 
only  in  the  equation  for  du/dt  and,  therefore,  leads  to  no  physical  paradoxes  in  the  limit  of 
circular  orbit.  However,  for  finite  values  of  eccentricity,  this  term  contributes  to  the  periapse 
precession. 

Another  seemingly  calamitous  thing  is  the  divergence  in  (82).  This  divergence,  however, 
entails  no  disastrous  physical  consequences,  because  the  term  (82)  shows  up  only  in  the 
planetary  equation  for  dM0/dt  and  simply  leads  to  a  steady  shift  of  the  initial  condition 
M0  . 


4.2  The  case  of  a  constant  precession  rate. 

The  situation  might  simplify  very  considerably  if  we  could  also  assume  that  the  precession 
rate  ft  stays  constant.  Then  in  equations  (69  -  74),  we  would  take  ft  out  of  the  angular 
brackets  and  proceed  with  averaging  the  expressions  (  (df/dCj)  x  g  —  /  x  ( dg/dCj )  ) 

only  (while  all  the  terms  with  ft  will  now  vanish).  It  is,  of  course,  well  known  that  this  is 
physically  wrong,  because  the  planetary  precession  has  a  continuous  spectrum  of  frequencies 
(some  of  which  are  commensurate  with  the  orbital  frequency  of  the  satellite).14  Nevertheless, 
for  the  sake  of  argument  let  us  go  on  with  this  assumption. 

Averaging  of  (75)  and  (80)  is  self-evident.  Averaging  of  (76  -  79)  is  lengthy  and  is  pre¬ 
sented  in  the  Appendix.  All  in  all,  we  get,  for  constant  ft : 


*  ■  <  I  fa  X  *  ~  *  *  fa  I  >  °  *  •  I  £  x  g  -  f  x  g  I  °  r  ■  ■  - 1  -  '  .  (87) 
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)  =  0  ,  Cj  =  e ,  ft ,  u ,  i ,  Ma 


(88) 


Since  the  orbital  averages  (88)  vanish,  then  e  will,  along  with  a  ,  stay  constant  for  as  long 
as  our  approximation  remains  valid.  Besides,  no  trace  of  ft  will  be  left  in  the  equations  for 
O  and  i  .  This  means  that,  in  the  assumed  approximation  and  under  the  extra  assumption 
of  constant  ft ,  the  afore  quoted  analysis  (23  -  36),  offered  by  Goldreich  (1965),  will  remain 
valid  at  not  time  scales  which  are  not  too  long.  At  longer  scales  (of  order  dozens  to  hundred 
of  millions  of  years)  one  has  to  take  into  consideration  the  back  reaction  of  the  short-period 
terms  upon  the  secular  ones  (Laskar  1990).  Beside  the  latter  issue,  the  problem  with  this 
approximation  is  that  it  ignores  both  the  long-term  evolution  of  the  spin  axis  and  the  short¬ 
term  nutations.  For  these  reasons,  this  approximation  will  not  be  extendable  to  long  periods 

14  The  case  of  the  Earth  rotation  and  precession  is  comprehensively  reviewed  by  Eubanks  (1993).  The 
Martian  short-time-scale  rotational  dynamics  is  of  an  equal  complexity,  even  though  Mars  lacks  oceans  and 
the  coupling  of  its  rotation  with  the  atmospheric  motions  is  weaker  than  in  the  case  of  the  Earth.  (Defraigne 
et  al  2003,  Van  Hoolst  et  al  2000,  Dehant  et  al  2000) 
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of  Mars’  evolution.  This  puts  forward  the  bigger  question:  what  maximal  amplitude  of 
obliquity  variations  could  Mars  afford,  to  keep  both  its  satellites  so  close  to  its  equatorial 
plane? 

Even  in  the  unphysical  case  of  constant  p  the  averaged  equation  (74)  for  the  osculating 
Ma  differs,  already  in  the  first  order  over  p  ,  from  equation  (48)  for  the  contact  Ma  . 
In  the  realistic  case  of  time-dependent  precession,  the  averages  of  terms  containing  p  do 
not  vanish  (except  p  ■  (  ( df/0Mo )  xg  —  /  x  (dg /dM0)  j  which  is  identically  nil).  These 
terms  show  up  in  all  equations  (except  in  that  for  a  )  and  influence  the  motion.  They  will 
be  the  key  to  our  understanding  the  long-term  satellite  dynamics,  including  the  secular  drift 
of  the  orbit  plane,  caused  by  the  precession  p. 


5  An  outline  of  a  more  accurate  analysis.  Resonances 
between  the  planetary  nutation  and  the  satellite  or¬ 
bital  frequency. 

Precession  of  any  planet  contains  in  itself  a  continuous  spectrum  of  circular  frequencies 
involved:15 


/OO 

P(s)  e~lst  ds  ,  (89) 

-OO 

some  modes  being  more  prominent  than  the  others.  For  our  present  purposes,  it  will  be 
advantageous  to  express  the  precession  rate  as  function  of  the  satellite’s  true  anomaly: 

/OO 

p(W)  e-lWu  dW  ,  (90) 

-OO 


IT  being  the  circular  “frequency”  related  to  the  true  anomaly  u  .  Needless  to  say, 
p(t)  p(v) ,  p(s) ,  and  p{W)  are  four  different  functions.  However,  we  take  the  lib¬ 
erty  of  using  the  same  notation  /7( . . . )  because  the  argument  will  always  reveal  which 
of  the  four  functions  we  imply.16  As  ever,  it  is  understood  that  the  physical  content  is 
attributed  to  the  real  parts  of  p(t)  and  p{v) . 

If  we  now  plug  the  real  part  of  (90)  into  (75  -  87)  and  carry  out  averaging  in  accordance 
with  formula  (109)  of  the  Appendix,  we  shall  see  that  the  secular  parts  of  these  /7-dependent 
terms  do  not  vanish.  They  reveal  the  influence  of  the  planetary  precession  upon  the  satellite 
orbital  motion.  Especially  interesting  are  the  resonant  contributions  provided  by  the  integer 
IT’s,  because  these  nutation  modes  are  commensurate  with  the  orbital  motion  of  the  satellite. 

For  plugging  p{y)  into  the  terms  (75  -  86),  it  will  be  convenient  to  rewrite  (90)  as 


POO 

PM  =  / 

Jo 


p{s\W)  sin(H'V)  +  pic\W)  cos(ITV) 


dW 


(91) 


15  A  more  honest  analysis  should  take  into  account  also  the  direct  dependence  of  the  planet’s  precession 
rate  upon  the  instantaneous  position(s)  of  its  satellite(s):  p,  =  /I(t;  a,  e,  ,  u> ,  i,  M0 )  .  However,  the 
back-reaction  of  the  satellites  upon  the  primary  is  known  to  be  an  effect  of  a  higher  order  of  smallness  (Laskar 
2004),  at  least  in  the  case  of  Mars;  and  therefore  we  shall  omit  this  circumstance  by  simply  assuming  that 

P  =  /7(f)  =  /7  ( t{v) ). 

16  Evidently,  p(v)  is  a  short  notation  for  /t(  t(v) )  .  It  is  also  possible  to  demonstrate  that,  in  the  limit 
of  vanishing  eccentricity,  ft(W)  ss  np(Wn)  and  W  ss  s/n  . 
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Insertion  of  the  latter  into  (75  -  86)  will  still  result  in  extremely  sophisticated  integrals.  For 
example,  the  term  (75)  will  have  the  following  orbital  average: 


(  P- 


3  iGm  (1  —  e2) 
2  V - a - 


(92) 


_3_  fGm  /  2y  f°°dw  f2w  du  l-i±\W)  sin(Wv)  +  ^(W)  cos(Wv) 
47 r  V  a  '  6  '  Jo  Jo  (1  +  e  cosv)2 


(the  averaging  rule  (109)  being  employed).  Evaluation  of  the  two  integrals  emerging  in  this 
expression  can,  in  principle,  be  carried  out  in  term  of  the  hypergeometric  functions,  but  the 
outcome  will  be  hard  to  work  with  and  hard  to  interpret  physically.  Even  worse  integrals 
will  show  up  in  the  averages  of  (76  -  86). 

To  get  an  idea  as  to  how  much  the  terms  (75  -  86)  contribute  to  the  secular  drift  of 
the  satellite  orbits  from  the  planet’s  equator,  it  may  be  good  to  first  calculate  these  term’s 
averages  under  the  assumption  of  small  eccentricity.  Then,  for  example,  (92)  will  simplify 
to 
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We  immediately  see  that  the  above  expression  consists  of  two  parts:  the  non-resonant  one 
and  the  resonant  one. 

This  topic  will  be  addressed  in  our  subsequent  paper  where  we  shall  try  to  determine 
how  much  time  is  needed  for  these  subtle  effects  to  accumulate,  to  provide  for  a  substantial 
secular  drift  of  the  orbit  plane. 


6  Conclusions 

In  this  article  we  have  prepared  an  analytical  launching  pad  for  the  research  of  long-term 
evolution  of  orbits  about  a  precessing  oblate  primary.  This  paper  is  the  first  in  a  series  and  is 
technical,  so  we  deliberately  avoided  making  whatever  quantitative  estimates,  leaving  those 
for  the  next  part  of  our  project. 

The  pivotal  question  emerging  in  the  context  of  this  research  is  whether  the  orbital  planes 
of  near-equatorial  satellites  will  drift  away  from  the  planetary  equator  in  the  cause  of  the 
planet’s  obliquity  changes.  Several  facts  have  been  established  in  this  regard. 

First,  the  planetary  equations  for  osculating  elements  of  the  satellite  do  contain  terms 
responsible  for  such  a  drift.  These  terms  contain  inputs  of  the  first  order  and  of  the  second 
order  in  p  ,  and  of  the  first  order  in  p ,  where  p  is  the  precession  rate  of  the  primary. 
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Second,  the  first-order  (but  not  the  second-order)  terms  average  out  in  the  case  of  a 
constant  precession  rate,  which  means  that  in  this  case  their  effect  will  accumulate  only  over 
extremely  long  time  scales.  We  would  remind  that  the  sort-period  terms  of  the  planetary 
equations  do  exert  back-reaction  upon  the  secular  ones.  While  in  the  artificial-satellite  sci¬ 
ence,  which  deals  with  short  interval  of  time,  the  short-period  terms  often  may  be  omitted,  in 
long-term  astronomical  computations  (dozens  of  million  years  and  higher),  the  accumulated 
influence  of  short-period  terms  must  be  taken  into  account.  A  simple  explanation  of  how 
this  should  be  done  is  offered  in  section  2  of  Laskar  (1990).  Besides,  the  contribution  of  the 
secular  second-order  terms  will,  too,  accumulate  over  very  large  time  intervals. 

Third,  the  first-order  drift  terms  do  not  average  out  in  the  case  of  variable  precession. 
Under  these  circumstances  they  become  secular.  This  means,  for  example,  that  the  turbulent 
history  of  Mars’  obliquity  -  history  which  includes  both  long-term  changes  (Ward  1973,  1974, 
1979;  Laskar  &  Robutel  1993;  Touma  &  Wisdom  1994)  and  short-term  nutations  (Dehant 
et  al  2000;  Van  Hoolst  et  al  2000;  Defraigne  et  al  2004)  -  might  have  lead  to  a  secular 
drift  of  the  initially  near-equatorial  satellites.  If  that  were  the  case,  then  the  current,  still 
near-equatorial,  location  of  Phobos  and  Deimos  may  lead  to  restrictions  upon  the  rate  and 
amplitude  of  the  Martian  obliquity  variations.  To  render  a  judgement  on  this  topic,  one 
should  compute  how  quickly  this  drift  accumulates.  Very  likely,  this  quest  will  demand 
heavy-duty  numerics. 
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Appendix. 

The  inertial  terms  emerging  in  the  planetary  equations. 


In  this  Appendix  we  write  down,  as  functions  of  the  true  anomaly,  the  ft-  and  ft- 
dependent  terms  in  (59  -  64).  We  also  calculate  their  orbital  averages  in  the  case  of  a 
constant  precession  rate  /7 . 


A.l.  The  basic  formulae 


Formulae  (55  -  60)  contain  the  two-body  unperturbed  expressions  for  the  position  and 
velocity  as  functions  of  the  time  and  the  six  Keplerian  elements,  (1)  and  (7).  To  find  their 
explicit  form,  one  can  employ  an  auxiliary  set  of  dextral  perifocal  coordinates  q  ,  with  an 
origin  at  the  gravitating  centre,  and  with  the  first  two  axes  located  in  the  plane  of  orbit: 


q  =  {  r  cos  v  ,  r  sin  u  ,  0  }  T 

The  corresponding  velocities  will  read: 


1  —  e2 
e  +  cost' 


{  cos  v  ,  sin  v  ,  0  }  T  .  (94) 


n  a  sin  v  n  a  (e  +  cos  v) 


VT 


n  a 
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{  —  sin  v  ,  (e  T  cos  v)  ,  0  } 


(95) 


the  radius  in  (94)  being  a  function  of  the  major  semiaxis  a  ,  the  eccentricity  e  ,  and  the 
true  anomaly  v  : 


1  -  e2 
1  +  e  cosu 


(96) 


the  true  anomaly  v  itself  being  a  function  of  a  ,  e  ,  and  of  the  mean  anomaly  M  = 
Ma  +  //o  n  dt  ,  where  n  =  (Gm)1/2  a-3/2  .  Then,  in  the  two-body  setting,  the  position 
and  velocity,  related  to  some  fiducial  inertial  frame,  will  appear  as: 


r  =  f  (Cj,...,Ce,  t)  =  R(fi,  i,  to)  q (a,  e,  Ma,t)  , 
r  =  ,g  (Ci, ...,  C6,  t )  =  R(0,  *,  u)  q(a,  e,  M0,t)  , 


(97) 


R(fi,  i,  u)  being  the  matrix  of  rotation  from  the  orbital-plane-related  axes  q  to  some  fixed 
Cartesian  axes  {x\ ,  x2,  Xs )  in  the  said  fiducial  inertial  frame  wherein  the  vectors  r  and 
r  are  defined.  The  rotation  is  parametrised  by  the  three  Euler  angles:  inclination,  i  ;  the 
longitude  of  the  node,  Q,  ;  and  the  argument  of  the  pericentre,  u  .  Thence,  as  well  known 
(see,  for  example,  Morbiclelli  (2002),  subsection  1.2), 
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cos  Cl  cos  lo  —  sin  12  sin  lo  cos  i  —  cos  12  sin  lo  —  sin  12  cos  lo  cos  i 

R  sin  12  cos  lo  +  cos  12  sin  lo  cos  %  —  sin  Cl  sin  lo  +  cos  Cl  cos  lo  cos  i 

sin  lo  sin  %  cos  lo  sin  i 

insertion  whereof,  together  with  (94),  into  the  first  equation  of  (97)  yields 
1  -  e2 

fi  =  a -  [  cos  12  cos  (u;  +  v)  —  sin  12  sin  (a;  +  u)  cos  i  ]  , 

1  +  e  costs 

1  -  e2 

f2  =  a -  [  sin  Cl  cos  (lo  +  is)  +  cos  fi  sin  (a;  +  is)  cos  i  }  , 

1  +  e  cos  is 

1  -  e2 

fs  =  a -  sin  (lo  +  is)  sin  i 

1  +  e  cosz^ 

Similarly,  substitution  of  (95)  and  (98)  into  the  second  equation  of  (97)  entails: 
ti  a 

gi  =  ,  [  —  cos  12  sin(cu  +  is)  —  sin  12  cos(u;  +  is)  cos  i  + 

VI  -  e2 


e  (  —  cos  Cl  sincu  —  sin  12  cos  a;  cosi)  ]  , 


n  a 


g2  =  7  [—  sin  12  sin  (a;  +  u)  +  cos  12  cos(uz  +  is)  cos  i  + 

Vl  -  e2 


e  (  —  sin  12  sin  cu  +  cos  12  cos  a;  cosi)  ]  , 


sin  Cl  sin  i 

—  cos  Cl  sin  i 

(98) 

cos  i 

(99) 

(100) 

(101) 


(102) 


(103) 


n  a  .  / 

gz  =  ,  smz  [  cos(tn  +  is)  +  e  costu 

Vl  -  e2 


(104) 


the  subscripts  1,2,3  denoting  the  x\ ,  x2  ,  Xz  components  in  the  fiducial  inertial  frame 
wherein  (97)  is  written. 


A. 2.  The  averaging  rule. 


From  equations 

cos  E 


it  follows  that 


e  +  cos  is 
1  +  e  cos is 


.  , ,  Vl  —  e2  sin  is 

sin  E  =  - 

1  +  e  cos  v 


dE  _  Vl  -  e2 
dv  1  +  e  cos  v 


(105) 


(106) 
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^From  the  first  of  formulae  (105)  and  from  the  Kepler  equation  one  can  derive: 


d  M  _  1  -  e2 

d E  1  +  e  coszz 

Together,  (106)  and  (107)  entail: 

dM  _  dM  dE  _  (1  -  e2)3/2 
du  dE  du  1  +  e  coszz 


whence 


1 

2~7T 


(1  -  e2f2 
2t r 


du 

(1  +  e  cos  u)2 


(107) 


(108) 


Calculation  of  the  integral  shows  that  the  right-hand  side  of  the  above  equation  is  equal 
to  unity,  which  means  that  the  secular  parts  should  be  calculated  through  the  following 
averaging  rule: 


(  •••  ) 


(1  -  e2)3/2  r 2” 
2  7T  JO 


du 

(1  +  e  cosz')2 


(109) 


In  what  follows  below,  we  try  hard  to  squeeze  the  calculations  as  much  as  possible,  but  at 
the  same  time  to  leave  them  verifiable  for  the  interested  reader. 


A. 3.  Calculation  of  p 


df  _  =*  dg 

Ta  X  g  -  /  X  Ya 


As  evident  from  (94)  and  from  the  first  equation  of  (97), 


df  df\  df\  / du\  f  df  dt  ( du\ 

da  \  da  J  \  du  )  \da ) .  ,,  a  dt  du  \  da 

\  J  v  \  /  0  ^  '  t,e,M0  \  /  t,e,  M0 


f  _  dt_  (dv' 
a  ^  ^  du  1  da  , 


t,  e,  M0 


and,  therefore, 


df  _  1 

da  a 


Similarly,  from  (95)  and  the  second  equation  of  (97)  it  ensues  that 


dg  =  (dg\  ( dg\  (dv' 

da  \  da  j  l  du  J  \  da  , 


t,  e,  M0 


(ill) 


(112) 


g  dg  dt  ( du  \ 

2  a  dt  du  1  da  J  , 

\  /  t,  e,  M0 


A  + 

2a 


Gm  dt  (du' 
|j|3  J  du  \(9ay 


t,  e,  M0 
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(110) 


wherefrom 


1 

2  a 


f  x  g  . 


(113) 


In  the  undisturbed  two-body  problem,  /  x  g  is  the  angular  momentum  (per  unit  of  the 
reduced  mass)  and  is  equal  to  J G  rn  a  ( 1  —  e2)  w? ,  where  the  unit  vector 


w  =  xq  sin  i  sin  —  x2  sin  i  cos  fi  +  x3  cos  i 


(114) 


is  normal  to  the  instantaneous  osculating  ellipse,  the  unit  vectors  xi,  x2,  X3  making  the 
basis  of  the  co-precessing  coordinate  system  xj ,  xg,  X3  (the  axes  xi  and  xg  belonging  to 
the  planet’s  equatorial  plane). 

Together,  (111)  and  (113)  will  give: 

df  -  «9g  3  =»  3  r- - — - —  _ 

—  xg-/x  —  =  — /  xg=  —  JGm  a  (1  -  e2)  w 
da  da  2  a  2  a  v 

(115) 


3  G  m  (1  —  e2) 


xi  sin  i  sin  fl  —  x2  sin  i  cos  Q  +  X3  cos  i 


and,  thereby, 


df  ,  dg  3  G  rn  (1  —  e2) 

"  '  1  Ta  X  g  -  1  X  d~a  =  2  ^  V - a - 


(116) 


where 


H±  =  jj,\  sin  i  sin  fl  —  fi2  sin  i  cos  Q  +  Hz  cos  i 


(117) 


Since,  for  constant  ft  ,  (116)  is  ^-independent,  then  in  the  uniform-precession  case  it  will 
coincide  with  its  orbital  average. 


A. 4.  Calculation  of  ft-\^-xg-fx^$ 

de  de 


Just  as  in  the  preceding  subsection,  one  can  write: 


df  =  df 

de  1  de 


df 

du 


2  e  +  cos  v  +  e2  cos  v 
(1  +  e  cos  v)  (1  —  e2) 


<9^ 

de ) 

j 

t,  a,  M0 

dt 

(du\ 

du 

Ue) 

df\  ,  (?a' 

de  8  du  [de  , 


t,  a.  M0 


(118) 


t,  a,  M0 
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whence 


df 


x  g 


2  e  +  cos  v  +  e2  cos  v 


/  x  g  ■ 


de  °  (1  +  e  cos  i/)  (1  —  e2) 

With  (95)  taken  into  account,  the  derivative  of  the  two-body  velocity  will  take  the  form: 


(119) 


dg  / dg  \  ( <>g\  (<>'■'' 

de  \  de  )  \dv)  \  <9e 


e  ,  dg  dt  ( dv\ 

=  - g  +  b  H - -  —  — 

M  1  —  e2  dt  dv  \de  j  M 

t j  a,  M0  \  /  t,  a,  M0 


e  ->  /  G  m  A  <9t  / <9zA 

—  g  +  b  +  -  —  /)  —  (  — ) 


1  —  e 


i/i3  ;  ^  v^eA,a, 


where 


(120) 


b  = 


(121) 


n  a 


VT^ 


{  —  cos  sin  u  —  sin  Cl  cos  a;  cos  %  ,  —  sin  sin  a;  +  cos  cos  a;  cos  z  ,  sin  i  cos  u  }  T  . 


From  this  we  easily  get,  with  aid  of  (114): 


/  x  b  = 


n  a2  y/1  —  e 2 
1  +  e  cos  v 


{  sin  Cl  sin  i  cos  v  ,  —  cos  Cl  sin  i  cos  v  ,  cos  z  cos  v  }  T  = 


(122) 


n  a ‘ 


VT 


1  +  e  cos  v 


cos  v  {  sin  i  cos  Q  ,  —  sin  i  cos  Cl  ,  cos  i  }  T  =  w 


n  a2  Vi  —  e2 
1  +  e  cos  v 


cos  Z'  . 


Together,  (119),  (120),  and  (122)  yield: 


df  ->  dg 

de  X  S  f  X  de 


2  e  +  cos  ^  +  e2  cos  v  e 

(1  +  e  cos  i/)  (1  —  e2)  1  —  e2 


/  x  g  -  /  x  b 


3  e  +  cos  v  +  2  e2  cos  rz - ~ ~  _  n  a2  a/1  —  e2 

—  m  - —  \/Gma  (1  —  e2)  —  w  — — : -  cos  v 


(1  +  e  cos  v)  (1  —  e2) 


1  +  e  cos  v 


—  w 


na 2  (3  e  +  2  cos  v  +  e2  cos  v) 


(1  +  e  cos  v)  a/1  —  e2 


whence 


(123) 


-  ,  df  _  -  dg 

/‘•  hr  x  g  -  /  x  ]r 

de  de 


A‘_l 


n  a2  (3  e  +  2  cos  z'  +  e2  cos  v) 
(1  +  e  cost')  \/l  —  e2 


With  aid  of  integrals 


r2ir 

Jo 


dv 


2  +  e2 


7T 


o  (1  +  e  cosz^)  (1  —  e2) 


,2V5/2 


(125) 
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r2?r  cos  is  dis 


-  3e 


(1  +  e  cos  is)3  (1  -  e2)'~ 


(126) 


it  is  easy  to  demonstrate  that,  under  the  assumption  of  constant  /2  ,  the  orbital  average  of 
(124)  is  nil. 


A. 5.  Calculation  of  ft  ■  (  x  g  -  /*  x  ^ 

\  Ouj  Ouj 


A  straightforward  calculation,  based  on  formulae  (94),  (95),  (97),  and  (98),  gives: 


OJ  2  Of3 

Sh  —  tt-  9  2 

OUJ  OUJ 


n  a 2  VI  —  e2 
1  +  e  costs 


e  sinil  sin*  sin^  ,  (127) 


df  -+\  df3  dfi  na2  V 1  -  e2  .  .  . 

^-xg  =  9i  ~e  9s  =  — - - e  cosil  sin!  sm  is  ,  (128) 

Ouj  J  Ouj  ouj  1  +  e  cos  is 


df  \  dfi 
X  g  = 

/  3 


dfz  na2  V 1  -  e2  .  . 

— —  gi  = - e  cos*  sm  is  .  (129) 

Ouj  1  +  e  cost' 


The  two-body  angular  momentum  cannot  depend  on  the  argument  of  the  pericentre, 


(/ x  g)  =  0  > 


(130) 


fx  f 


(131) 


whereby 


s  x  M 

du 


'll  Cl2  \f\  —  P2 

—  2  - -  e  sin Q  sin*  sinz/  ,  (132) 

1  +  e  cos  is 


df  =*  dg 

dUxe~fx-K, 


'll  (l2  \/l  —  p2 

2  - e  cos  sin  *  sin  is  , 

1  +  e  cos  is 


(133) 


df  ^  ?  dg 

a;  x  s  -  /  x  K, 


n  a 2  a/1  —  e2 

—  2  - e  cos  *  sm  z/ 

1  +  e  cos  z^ 


(134) 
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Altogether,  these  formulae  may  be  written  down  as: 


df  =*  5g 

to  x  6  “  f  x  to 


—  2  w 


n  a 2  a/1  —  e2 
1  +  e  cost' 


e  sin  v 


(135) 


w  being  the  unit  vector  pointing  in  the  direction  of  the  angular  momentum.  It  is  given  by 
(114).  Finally, 


_  .  df  _  -*  5g 

1  to  x  g  “  f  x  to 


—  2 


n  a 


2  \J  1  -  *2 


1  +  e  cos  v 


e  sin  v 


(136) 


Since  sint'  is  odd  and  (1  +  e  cosu)  3  is  even,  it  is  obvious  that,  for  a  constant  p  ,  the 
orbital  average  of  (136)  will  vanish. 


A. 6.  Calculation  of  p 


df 

512 


X  g 


Once  again,  formulae  (94),  (95),  (97),  and  (98)  yield: 


df 


df, 


dfs 


Qn  ^  8  qo  on  92 


512 


n  a‘ 


512 


512 


(137) 


^T 


1  +  e  cosu 


cosl 2  cos(o;  +  u)  —  sin  12  sin(<u  +  v)  cos  i }  cos  (a;  +  u)  sin  i 


+ 


n  a " 


VT 


1  +  e  cos  nu 


e  [cos  £2  cos(cu  +  v)  —  sin  12  sin(a;  +  u)  cos  i  ]  cosu;  sin  i  , 


df 

512 


X  g 


dfs  dfi 

512  91  512  93 


(138) 


n  a2  \/l  —  e2 
1  +  e  cos  v 


[sin  12  cos(cu  -|-  v)  |  cos  12  sin(cu  +  v)  cos  i  ]  cos  (a;  +  v)  sin  i 


+ 


n  a 2  \/l  —  e2 
1  +  e  cos  v 


e  [sin  12  cos(tu  +  u)  +  cos  12  sin(a;  +  v)  cos*]  cos  a;  sin*  , 


df  A  dh  df2 

512  X  S  512  92  512  91 


(139) 
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na2\/l  —  e2 
1  +  e  cos  v 


{  [  —  sin  ft  cos(o;+  v) 


—  cos  ft  sin(u;  +  v)  cos  i  ] 


[  —  sin  ft  sin  (a;  +  v)  +  cos  ft  cos  i  cos(u;  +  v)\  — 


[cosfl  cos((n  +  z/  —  sinfl  sin(a;  +  v)  cosi]  [—  cosQ  sin  (a;  +  v)  —  sin  11  cos  i  cos(a;  +  v)\  } 


n  a 2  a/1  —  e2 

+  — -  e 

1  +  e  cos  v 


{[ 


sinl!  cos(tn  +  v)  — 


cos  11  sin(a;  +  v)  cos  i  )  [—  sin  11  sin  a;  +  cosH  cos  i  cos  a;] 


—  cosH  cos(a;  +  u) 


sin  H  sin(u;  +  v)  cos  i 


—  cos  H  sin  a;  —  sin  H  cos  i  cos  u  ]  }  = 


na2V  1  —  e2 
1  +  e  cos  v 


|  sin(cu  +  u)  cos(lv  +  u)  sin2  i  +  e  sinw  cos(ui  +  u)  —  sin(a;  +  v)  cos  a;  cos2  %  j  , 


ag  og3  a  g2 

f  sn  ,  an  k  on 


(140) 


n  a2  \/l  —  e2 
1  +  e  cos  v 


[—  sin(cn  +  v)  sin  i  )  [—  cos  ft  sin(a;  +  v)  —  sin  11  cos(u;  +  u)  cos*] 


+ 


n  a 


vT~ - 


1  +  e  cos v 


e  [  —  sin  (a?  +  v)  sin  i  ]  [  —  cos  ft  sin  u>  —  sin  ft  cos  uj  cos  i  ] 


n  a 2  y/l  —  e2 
1  +  e  cost' 


sin  % 


sin2  {uj  +  v)  cos  ft  +  sin  (u  +  u)  cos  (u  +  v)  sin  ft  cos  i 


+ 


n  a2  a/I  —  e2 
1  +  e  cos  nu 


e  sin  i  [  sin  (to  +  v)  sin  u  cos  ft  +  sin  (u ;  +  v)  cos  uj  sin  ft  cos  i  ]  , 
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( 


<9g 

99 


9gi  9g3 
73  <90  71  99 


(141) 


n  a2  a/1  —  e2 
1  +  e  coszv 


sin(u;  +  z/  sin  *  [sin  9  sin(u;  +  v) 


cos  9  cos(u;  +  v)  cos  *  ] 


fl 

+  -  e  sin^  +  z/  sin  *  [sin  9  sin  a; 

1  +  e  cos  v 


cos  9  cos  lo  cos  i  ] 


( 


df 

99 


x 


dhn 

99  92 


no?\J\  —  e2 
1  +  e  cos  v 


{  [  cos  9  cos(a;  +  v) 


(142) 


—  sin  9  sin(c<;  +  v)  cos  i  ]  [  —  cos  9  sin  ( lo  +  u) 


sin  9  cos*  cos(u;  +  ^)]  — 


sin  9  cos(a;  +  v)  +  cos  9  sin(a;  +  u)  cos  *  ]  [  sin  9  sin  (lo  +  u)  —  cos  9 


COS  *  cos 


[u  +  I/)]  } 


n  a2  a/1  -  e2  ,  r  _  , 

H - e  {  [cos 9  cos(a;  +  v) 

1  +  e  cos  v 


sin  9  sin  (a;  +  v)  cos  *  ]  [—  cos  9  sin  a;  —  sin  9 


cos  *  cos  UJ 


—  [sin  9  cos(t<;  +  v)  +  cos  9  sin(a;  +  v)  cos  *  ]  [sin  9  sin  a;  —  cos  9  cos*  cos  a;]  } 


na2  a/1  —  e2 
1  +  e  cos  v 


|  —  sin(a;  +  v)  cos(u;  +  z/)  sin2*  +  e  —  sin  lo  cos(u;  +  v)  +  sin(u;  +  v)  cosu;  cos2* 


In  (137  -  138)  and  in  (140  -  141)  we  benefitted  from  f3  being  independent  of  9  .  Now,  by 
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subtracting,  accordingly,  (140)  from  (137),  and  (141)  from  (138),  and  (142)  from  (139),  we 
arrive,  after  some  intermediate  algebra,  at  the  following  three  expressions: 


9f  ?  <9g 

an  x  s  f  x  on 


(143) 


n  a2  a/1  —  e2 
1  +  e  cost' 


sin  i  {  cos  Q  cos  [  2  (uj  +  u)  ]  —  sin  Q  cos  i  sin  [  2  (a;  +  u)  ]  } 


+ 


n  a2  a/1  —  e2 
1  +  e  cos  v 


e  sin*  {  cosfl  cos  (u  +  2  a>)  —  2  sinQ  cos  i  sin  (uj  +  u)  cos  a;  }  , 


df  ?  dg 

an  x  s  f  x  an 


(144) 


n  a ‘ 


1  +  e  cost' 


sin  i  {  sin  £2  cos  [  2  (iv  +  v) }  +  cos  £2  cos  i  sin  [  2  (u  +  u)  ]  } 


Tl  Qj2  a/T  — 

+  - - - - —  e  sini  {  sin  £2  cos  (u  +  2  iv)  +  2  cos  £2  cos  *  sin  (u  +  v)  cos  a;  }  , 

1  +  e  cos  u 


df  _  -  <9g 

mxe~f)<xi 


(145) 


n  a 


1  +  e  cosr/ 


sin2  i  sin  [  2  (u  +  z')  ] 


n  a2  a/1  —  e2 
1  +  e  cost' 


e  |  2  sin  cu  cos(u  +  u)  —  2  sin(a;+z/)  cos  a;  cos2  i  j 


nar  a/1  —  e 


1  +  e  cos  i/ 


sin2  %  sin  [  2  (tu  +  z/)  ]  + 


tz(z2a/1  —  e2 
1  +  e  cos  z/ 


2  e  |  —  sin  v  +  sin2  i  cos  u  sin  (uj  +  v)  } 


Hence,  the  overall  expression  for  ft  ■  ( ( df/dfl )  x  g  —  /  x  (<9g/<9£2))  will  be  given  by 
(78).  To  average  this  expression,  in  the  simple  case  of  a  constant  ft  ,  it  will  be  sufficient 
to  separately  average  the  above  three  expressions  (143  -  145)  and  then  to  multiply  the 
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averages  by  /;,]  ,  fi2  and  //3  ,  correspondingly,  and  to  add  everything  up.  To  carry  out  this 
procedure,  it  will  be  convenient  first  to  regroup  terms  in  (143  -  145)  and  to  through  out  the 
inputs  proportional  to  the  odd  functions  sin  v  and  sin  2v  ,  because  these  inputs  will  vanish 
after  averaging  via  formula  (109).  Thus  we  get: 


df  _  - 


X 


(146) 


n  a2  a/1  —  e 
1  +  e  cosz' 


2 


sin  i  {  cos  fi  cos  2a;  cos  2v 


cosfZ  sin  2a;  sin2z2 


—  sin  Q  cos  *  sin  2a;  cos  2v  +  sin  fZ  cos  *  cos  2a;  sin  2 v  } 


Tl  (X 2  \f\  ~  g2 

+  -  e  sin  *  {  cos  fZ  cos  2a;  cos  v 

1  +  e  cost; 


cos  f Z  sin  2a;  sin  u 


—  sin  fZ  cos  i  sin  2a;  cos  v  +  2  sin  fZ  cos  *  cos2  u>  sin  v 


.  n  a2  a/  1  —  e2  .  . 

(  -  sin  * 

1  +  e  cos  u 


(  cos  fZ  cos  2a;  —  sin  f Z  cos  i  sin  2a;  )  (  cos  2v  +  e  cos  v  )  )  = 


n  _  p2l3/2 

------ vr 


— k  .  .  ,  _  _  ^  ,  r2n  cos  2z/  +  e  cos  z/  , 

sini  cosiZ  cos  2a;  —  smSZ  cos*  sin  2a;  )  /  - n- dv  =  0  , 

2  7r  v  Wo  (1  +  ecosz/)3 


<9/  -  r  <9g  - 

on  x  s  ^  x  an  )  ^ 

2 


(147) 


n  a2  y/l  —  e2 
1  +  e  cos  v 


sin  i  {  sin  fZ  cos  2a;  cos  2 v  —  sin  fZ  sin  2a;  sin  2v 


+  cosfZ  cos*  sin  2a;  cos2z/  +  cosfZ  cos*  cos  2a;  sin  2v  } 
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Tl  (X2  ^/l  ~  g-2 

+  -  e  sin  i  {  sin  n  cos  2u  cos  u 

1  +  e  cos  v 


sin  n  sin  2a ;  sin  v 


+  cos  Cl  cos  i  sin  2a;  cos  v  +  2  cos  n  cos  i  cos2  a;  sin  v 


n  a 2  a/1  —  e2 
1  +  e  cos  zv 


sin  i  (sinO  cos  2a;  +  cosJl  cos  i  sin  2a;)  (cos  2 u  +  e  cosz/  )  = 


(1  -  e2)3/2  2 

2  7T 


e2  sin  i  (  sin  0  cos  2a;  +  cos  Cl  cos  i  sin  2a;  )  / 

ao 


2,r  cos  2z^  +  e  cos 

3~  dv  =  0  , 


(1  +  e  cos  z/)" 


<9/  r  |  v 

an  x  s  x  an 


(148) 


na2  a/1  —  e2 
1  +  e  cos  i/ 


sin2  i  (  sin  2a;  cos  2v  4-  cos  2a;  sin  2v  )  + 


na2  a/1  —  e2 
1  +  e  cos  i/ 


2  e  (  —  sin  z/  +  sin2  i  cos  a;  sin  a ;  cos  z/  +  sin2  i  cos2  a;  sin  z^  )  = 


( 


na2  a/1  —  e2 
1  +  e  cos  z/ 


sin2  i  sin  2a;  (cos  2z/  +  e  cos  u)  )  = 


n  _  g2)3/2  _ 

- - - n  a2  a/1  —  g2  sin2  i  sin  2a; 

2  7T 


cos  2z/  +  e  cos  z'  , 

— - -o-  dv 

(1  +  e  cos  v) 


0  . 


We  see  that  the  averages  of  all  three  Cartesian  components  of  (  (df  /dCl)  x  g  -  f  x  (dg/dCl) ) 
vanish.  In  the  developments  (146  -  148)  the  following  integrals  were  used: 


and 


cos  v 

(1  +  e  cosz/)3 


1 1  —  e 
1  +  e  (1 


3  7 re 

e)3  (1  +  e)2 


cos2z^ 

(1  +  e  cos  vf 


1  —  e  3  7r  e2 
1  +  e  (1  -  e)3  (1  +  g)2 


(149) 


(150) 
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A. 7.  Calculation  of  ft 


df  r  <9g 

Ttt  x  g  -  /  x  — 
Ol  Ol 


Just  like  the  preceding  subsections,  this  one  is  based  on  formulae  (94),  (95),  (97),  and 
(98): 


( 


df 

di 


x 


(151) 


no2  a/  1  —  e2 
1  +  e  cos  v 


sin(u;  +  u)  sin  i  cos  O  ]  [  sin  i  cos  (iv  +  u)  ]  — 


cos  i  sin(u;+i/)]  [—  sinfi  sin  (a;  +  v)  +  cosfi  cos(u;  +  v)  cos?]  } 


+ 


n  a2  a/1  —  e2 
1  +  e  cosu 


e  {[ 


sin(cj  +  v)  sin*  cos  14  ]  [sin?  cos  a;  — 


cos?  sin(a;  +  u)  ]  [—  sinSJ  sin  a;  +  cosU  cos  a;  cos?]  }  = 


na2\/  1  —  e2 
1  +  e  cos  u 


{ 


sin  fl  sin2  (u  +  v)  cos  i 


cos  Q  sin  (a;  +  v)  cos  (u>  +  v) 


+ 


e  sin  (u  +  v)  [  sin  O  sin  tu  cos  ?  —  cos  Q  cos  u)  ]  }  , 


( 


df 

di 


x 


(152) 


t?ci2a/1  —  e2 
1  +  e  cos  z/ 


sin(u;  +  ?/)  cos  i 


cos  Q  sin  (iv  +  v)  —  sin  Vt  cos  ?  cos  (w  +  v)  ]  — 
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[  sin  ft  sin  i  sin(a;  +  v)  ]  |  sin  i  cos  (10  f  v)  \  } 


n  a2  \/l  —  e2  ,  .  .  , 

H - e  {  [  sm(w  +  v)  cos  % 

1  +  e  cos  nu 


cos  ft  sin  lo  —  sin  ft  cos  lo  cos  i  }  — 


[  sin  ft  sin(a;  +  u)  sin  i  } 


[  sin  i  cos  lo  }  }  = 


nc?\J  1  —  e2 
1  +  e  cos  v 


sinfi  sinfw  +  v)  cos(a;  +  v) 


cos  ft  sin2  (a;  +  u)  cos  i 


+ 


e  [  —  sinfi  sin  (lo  +  v)  cos  a;  —  cosfi  cos*  sinu;  sin  (a;  +  u)  ]  }  , 


( 


df 

di 


x 


(153) 


m2V  1  —  e2 
1  +  e  cos  v 


{  [  sin(cu  +  u) 


sin  i  sin  ft 


sin  sin  (lo  +  u) 


+  cos  cos  i  cos(u;  +  v)  ]  — 


—  cos  ft  sin(c<;  +  u)  sin  i 


cos  ft  sin  (lo  +  u)  —  sin  ft  cos(a;  +  v)  cos  i  }  } 


+ 


n  a2  \/l  —  e2 
1  +  e  cos  v 


e 


{  [  sin(k;  4-  v)  sin  i  sin  ft 


sin  ft  sin  lo  +  cos  ft  cos  i  cos  lo  ]  — 


[—  cos  ft  sin(u;  +  i/)  sin  i  ]  [—  cosfl  sin  a;  —  sinfi  cos  a;  cos  i  ]  }  = 


na2\/ 1  —  e2 
1  +  e  cos  v 


|  —  sin2(a;  +  ^)  sin  i  —  e  sin(cu  +  z/)  sin  a;  sin*  |  , 


( 


<9g 

di 


f  d93  f  dg2 

J  2  r\  •  J  3  r\  • 

oi  a% 


(154) 
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n  a2  a/1  —  e2 
1  +  e  cos  v 


{  [sinfi  cos {uj  +  v)  +  cosQ  sin(a;  +  u)  cos?] 


cos(a;  +  v)  cos  i  — 


sin  (a;  +  u)  sin  i  [  —  cos  Q  cos {uj  +  u)  sin  i  ]  } 


+ 


n  a2  a/1  —  e2 
1  +  e  cos  v 


e  {  [sinfl  cos(a;  +  u)  +  cosfi  sin(u;  +  v)  cos?]  cos  u>  cos  i  — 


sin(a;  +  u)  sin  ?  [—  cosf2  cos  a;  sin  i  ]  }  = 


n  a 


1  +  e  cosu 


|  cos  SI  sin(u;H-i/)  cos(u;  +  v)  +  sin  12  cos2(a;  +  ^)  cos  ? 


+ 


e  [cos  SI  sin(u;  +  v)  cos  uj  +  sin  11  cos(w  +  v)  cos  a;  cos  i]  }  , 


( 


dg 

di 


dgi  dg3 

73  di  71  di 


(155) 


n  a2  a/1  —  e2 
1  +  e  cos  v 


{  sin  (a;  +  v)  sin  i 


sin  S4  sin  i  cos (u  +  v)  — 


[  cos  Vt  cos(a;  +  v)  —  sin  f2  sin(a;  +  u)  cos  i\  cos(u;  +  u)  cos  i  } 


'll  CL2  ^/l — C2 

H - e  {  sin(a;  +  z/  sin  i  sinfi  sin?  cos  a;  — 

1  +  e  cost' 


[cosS~2  cos(n;  +  v)  —  sinfi  sin(a;  +  v)  cos?]  cos  a;  cos?  }  = 


n  a2  a/I  —  e2 
1  +  e  cos  v 


sin  Q,  sin  (a;  +  u)  cos  (a;  +  v) 


cos  SI  cos2  (a;  +  v)  cos  ? 


+ 
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e  [sin  SI  sin  (a;  +  u)  cos  a;  —  cos  SI  cos(a;  +  u)  cos  a;  cos  i  ]  }  , 


f  x  5 i\  =  f  _  f.  I 

^  O'J  7 1  o  •  72  ^  • 

cn  /  3  m  a* 


(156) 


na 2  v  1  —  e2 
1  +  e  cos  ^ 


{  [  cos  SI  cos(u;  +  v)  —  sin  Cl  sin(u;  +  u)  cos  i  ]  [  —  cos  Cl  sin  i  cos(u;  +  v)  }  — 


[sin  SI  cos  (a;  +  v)  +  cos  SI  sin(a;  +  u)  cos  i\  sin  SI  cos  (u>  +  v)  cos  i  } 


na 2  \/l  —  e2 


H - e  {  [cos  SI  cos  (u>  +  u)  —  sin  SI  sin(^  +  z/)  cos  i  } 

1  +  e  cos  v 


cos  S2  sin  i  cos  u 


sin  SI  cos(a;  +  v)  +  cos  SI  sin  (a;  +  u)  cos  i]  sinS2  cos  a;  cos  i  }  = 


n  a2  VI  —  e2 
1  +  e  cosn 


|  —  sini  cos2(a;  +  v)  —  e  sin  i  cos(w  +  u)  cos  a;  } 


Now,  by  putting  together,  accordingly,  (151)  with  (154),  (152)  with  (155),  and  (153)  with 
(156),  we  arrive,  after  some  intermediate  algebra,  to  the  following  expressions  for  the  three 
Cartesian  components: 


df  ?  dg 

TV  x  6  -  f  x  7T 
oi  oi 


(157) 


na2\J  1  —  e2 
1  +  e  cos  v 


|  sin  Cl  cos  i  (  sin2(tu  +  v)  —  cos2(u  +  v)  )  —  2  cos  SI  sin  {uj  +  u)  cos  (u>  +  v)  + 


e  [  sin  Cl  cos  i  (  —  cos  (tu  +  v)  cos  u  +  sin  (u  +  v)  sin  uj  )  —  2  cos  Cl  sin  (w  +  v)  cos  uj  }  } 


na2  vl  —  e2 
1  +  e  cos  v 


|  (  —  sin  Cl  cos  i  cos  2u  —  cos  SI  sin  2u  )  (  cos2  u  —  sin2  v  j  -|- 
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2  sin  v  cos  v  (  sin  Cl  cos  i  sin  2o>  —  cos  Cl  cos  2o>  )  ]  + 


(  —  sin  Cl  cos  i  cos  2 o>  —  cos  Cl  sin  2 u  )  cos  u  +  ^  sin 


cos  i  sin  2a;  —  2  cos  Cl  cos2  a; )  sin  v 


^  sin^  |  , 


df  ?  dg 

-W7  X  g  -  /  x  — 
Ol  Ol 


(158) 


na2  a/1  —  e2 
1  +  e  cos  v 


|  cosfl  cosi  (  cos2(a;  +  u)  —  sin2(a;  +  ^)^—  2  sinQ  sin(a;  +  ^)  cos  (a;  +  v) 


+ 


e  [  cos  Cl  cos  i  (  cos  (a;  +  u)  cos  oj  —  sin  (a;  +  u)  sin  iv  )  —  2  sin  Q  sin  (a;  +  u)  cos  u  ]  } 


na2\/l  —  e2 
1  +  e  cos  v 


(  cos  Cl  cos  i  cos  2a;  —  sin  Cl  sin  2a; ) 


cos2  v 


sin2  v 


2  sin  v  cos  v  (  —  cos  (1  cos  i  sin  2a;  —  sin  cos  2a;  )  ]  + 


e  (  cos  Cl  cos  i  cos  2a;  —  sin  Cl  sin  2a; )  cos  v  +  (  —  2  sin  Cl  cos2  uj  +  cos  Cl  sin  2a;  cos  i  J  sin  v  |  , 


df  ?  dg 

77  X  g  -  /  x  — 
Ol  Ol 


(159) 


na2  yl  —  e2 


1  +  e  cos  v 


|sinz  cos2(a;  +  ^)  —  sin2(a;  +  ^) 


+ 


e  [sin*  cos  a;  cos(a;  +  v)  —  sin?  sin  a;  sin(a;  +  u)  ]  }  = 
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nc?\J  1  —  e2 
1  +  e  cos  v 


Sill  % 


cos 


2a;  ( 


cos2  v  —  sin2  v 


2  sin  2 uj  sin  is  cos  is 


+ 


e  [  sin  i  cos  2c o  cos  is  —  sin  i  sin  2a;  sin  v  ]  }  . 
Averaging  of  the  afore  calculated  three  Cartesian  components  of  l(df/di)  x  g  —  /  x  (<9g/<9 i 

almost  exactly  coincides  with  averaging  of  (  (df  /dQ,)  x  g  -  /x  (<9g/<9f2)^  presented  in 
the  preceding  subsection,  and  the  result  is  the  same,  nil.  Indeed,  as  can  be  easily  seen  from 
(157  -  159),  after  we  through  out  the  (vanishing  after  averaging)  inputs  proportional  to  the 
odd  functions  sin  v  and  sin  2^  ,  each  of  the  three  expressions  (157  -  159)  will  be  propor¬ 
tional  to  (cos2i/  +  e  cos^)/(l  +  e  cosv)  .  Averaging  of  this  expression  via  formula  (109) 
will  vanish,  as  can  be  easily  seen  from  (149  -  150). 


A. 8.  Calculation  of  /2 


df 

0Mo 


x  g  -  /  x 


dM0 


and 


This  calculation  will  be  the  easiest:  as 

df  =  df  dt 
dM0  ~  dt  dM0 

dg  dg  dt 


=  g 


dM0  dt  dM„ 


dt 

dMo 
G  m 

"  T7F 


/ 


dt 

dK 


then,  obviously, 


df 

dM0 


x  g  -  /  x 


dMa 


/x  -  (0  -  0)  =  0 


(160) 


(161) 


(162) 


A. 9.  Calculation  of  fi  -  f  x 


d_l 

da 


With  aid  of  (110)  and  (68)  we  establish: 
f  df  r  .  (  dt  \  (dE\ 

-  /x  Ta=  -/xg  (ssj  „  (srj,  M 

\  /  a,e,Mo  \  /  t,e,Mo 


(163) 


w  JCma  (1  —  e2)  - 

V  -71 


1  1  -  e2 


'dE' 


n  1  +  e  cos  is  \  da 


—  w  a 


t,  e,  M0 


2  (1  —  e2)3/2  (dE" 
1  +  e  cos  is  \  da  , 


t.  e.  M0 
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where  we  used  the  formula 


dt 

dE 


1  —  e  cos  E  1  1  —  e2 


a,  e,  M0 


n 


n  1  +  e  cos  v 


(164) 


that  follows  from  the  Kepler  equation  and  from  the  first  equation  of  (105). 
It  remains  to  compute  dE / da  .  From  the  Kepler  equation 


it  follows  that 


wherefrom 


E  —  e  sin  E  =  M0  +  n  (  t  —  tQ  ) 


dE  ( 1  —  e  cos  E)  =  —  -  n  a  1  (t  —  ta)  da 


(165) 


(166) 


'dE\ 
da  ) 

'  '  t,e, 


M0 


3  _j  n  (t  -  t0) 
2  1  —  e  cos E 


3  _1  E  —  e  smE  —  M0 

-  a  - 

2  1  —  e  cos E 


(167) 


3  _!  1  +  cos;/ 

n  a  — - —  [t  -  t0)  . 


2  1  -  e2 

Insertion  of  the  above  into  (163)  will  give: 

;r  df  3  _  ..  ,  .  n — 

-  f  x  —  =  -  w  an  (t  -  to)  vi  - 


da 


expression  linear  in  t  —  tQ  . 


A. 10.  Calculation  of  /!•  I  -  /  x  ^ 

de 


Repeating  the  line  of  reasoning  presented  in  the  preceding  calculation,  we  derive  from 
(118),  (63),  and  (164)  the  following: 


-*  df  ->  _  (  dt' 

'  /  X  =  “/Xg  (dE. 


‘dE\ 


—  wo 


\  de  I 

a,  e,  M0  \  /  t,  e,  Ma 

To  compute  dE /de  ,  one  can  start  with  the  Kepler  equation, 

E  —  e  sin  E  =  M0  +  n  (l  —  ta) 
which  yields,  for  t,  a,  Ma  being  fixed: 

dE  —  e  cos  E  dE  —  sin  E  de  =  0 


2  (1  ~  e2)3/2  fdE' 
1  +  e  cos  v  \  de 


t,  a,  M0 


(169) 


(170) 


(168) 
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insertion  of  (105)  wherein  will  entail 


All  in  all, 


d  E  \  sin  E  sin  v 

9 eJt,a,Mo  1  -  e  cos  E  y/l  -  e2 


—  w  a2 


(1  ~  e2) 

1  +  e  cos  v 


sin  v 


(171) 


(172) 


A. 11  Calculation  of  pt 


Evidently,  this  vector  is  perpendicular  to  the  instantaneous  orbital  plane  and  is  u- 
independent.  A  direct  computation,  based  upon  (99  -  101),  will  lead  us  to: 


9±  x  f 

du 


9h  ,  9h  , 

du  h  du  h 


(173) 


r2  [—  sin  12  sin(u;  +  v)  +  cos  12  cos(u;  +  v)  cos  *  ]  sin(o;  +  u)  sin  *  — 


r2  cos(u;  +  v)  sin  i  [sin  12  cos  (u  +  u)  -1  cos  12  sin(a;  +  i/)  cos  *  = 


2  (1  -  e2)2  r  .  . 

a  - ■ — -77  —  sin  12  sin  * 


(1  +  e  cos  vy 


9f 


eh  ,  _8h 

du  du  h 


(174) 


r2  cos  (a;  +  v)  sin  *  [  cos  12  cos(tu  +  v)  —  sin  12  sin(a;  +  v)  cos  i  ]  — 

r2  [—  cos  12  sin(u;  +  u)  —  sin  12  cos(cu  +  u)  cos*]  sin(a;  +  v)  sin*  = 

2  (1  -  e2)2  .  .  1 

a  - o  cos  1 1  sin  *  , 

(1  +  e  cos  u  f  1  J 


df 

du 


x 


dfi  ,  _  9_h 

du  2  du 


(175) 
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r2  [—  cosH  sin(u;  +  u)  —  sin  12  cos(uj  +  v)  cos  i  }  [sin  12  cos (oo  +  u)  +  cos  12  sin(a;  +  v)  cos  * 
r2  [—  sin  12  sin  (a;  +  v)  +  cos  12  cos(u;  +  z/)  cos  i  ]  [cos  12  cos(u;  +  v)  —  sin  12  sin(u;  +  u)  cos 


Briefly, 


2^2 


—  r2  cos  i  =  —  a2 


(1  -  e2) 


(1  +  e  cos  v)z 


cos  * 


di 

duj 


2\2 


/  ^  =  —  w  a2 


(1  -  e2) 


(1  +  e  cos  v)z 


(176) 


A. 12.  Calculation  of  jl  -/x 


df 

<912 


df  x  A  =  dh  f  _  dh  f  _ 
<912  f  on  '3  on  '2  ' 


r2  [cos  12  cos (u  +  v)  —  sin  12  sin(u;  +  z/)  cos  *  \  sin(a;  +  v)  sin  *  = 


(177) 


(1  -  e2) 


2^ 


(1  +  e  cos  v)‘ 


cos  12  cos  (cv  +  v)  —  sin  12  sin(u;  +  v)  cos  *  ]  sin(u;  +  u)  sin*  , 


df  ? 

an  x  f 


dkf_dhf 
on  Jl  on  13 


(178) 


—  r2  [—  sin  12  cos(c<;  +  v)  —  cos 12  sin(c<;  +  z/)  cos*]  sin(a;  +  v)  sin* 


(1  -  e2) 


2n2 


(1  +  e  cos  uf 


[  sin  12  cos(o;  +  v)  +  cos  12  sin(u;  +  v)  cos  *  ]  sin(u;  +  u)  sin  *  , 


?L  x  f  1  =  dfi  ,  _  d_h  , 

an  f  j  on h  on h 

/  3 


(179) 
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r2  [—  sin  SI  cos(u)  +  v)  —  cos  SI  sin(cj  +  u)  cos  i  ]  [sin  SI  cos(a;+^)  +  cos  SI  sin(a;  +  u)  cos  i 
r2  [  cos  Cl  cos(c(;  +  v)  —  sin  Cl  sin  (a;  +  u)  cos  i  ]  [  cos  Cl  cos(a;  +  u)  —  sin  Cl  sin  (a;  +  v)  cos  i  ] 


—  cos2(c^  +  ^)  —  sin2  (cj  +  i/)  cos2  % 


2\2 


—  a 


.2  (1  -  **) 


(1  +  e  cos  v)' 


cos2  (a;  +  v)  +  sin2  (a;  +  v)  cos2 


A. 13.  Calculation  of  pi  I  -fx 

di 


dJif_dA 

di  h  di 


(180) 


r 2  {  [  —  cos  Cl  sin(cj  +  v)  sin  i  ]  sin(oj  +  v)  sin  i  — 

sin(cu  +  v)  cos  i  [sin  SI  cos(a >  +  v)  +  cos  SI  sin(u;  +  v)  cos  i  ]  }  = 

(1  -  e2)2 

a2  - ~2  [  —  cos  SI  sin(u;  +  v)  —  sin  SI  cos(a;  +  u)  cos  i  ]  sin(u;  +  u)  , 

(1  +  e  cosu ) 


df 

di 


dh 

di 


(181) 


r2  {  sin(c<;  +  u)  cos  i  [  cos  Cl  cos(u  +  v)  —  sin  SI  sin(o;  +  v)  cos  i  ]  — 


sin  SI  sin(u;  +  v)  sin  i  sin(o;  +  v)  sin  i  }  = 


a 


2 


(1  ~  e2)2 
(1  +  e  cos  is)2 


sin  SI  sin(u;+z/)  +  cos  SI  cos(u;  +  v)  cos*]  sin(u;  +  v)  , 
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3 


d±  X  f 

di 


dfi  ,  df2 

~wr  J  2  ~  ~W7~  J 1 
Ol  Ol 


r2  [  sin  fl  sin(a;  +  v)  sin  i  ]  [  sin  Q  cos((u  +  v)  +  cos  fl  sin(o;  +  u)  cos  i  ]  — 


r2  [—  cos  12  sin(a;  +  ^)  sin  i  ]  [cos  12  cos(u>  +  i2)  —  sin  12  sin(a;  +  v)  cos  i  ] 


(182) 


r2  sin(u;  +  u)  cos(a;  +  v)  sin  i  =  a2 


(1  -  e2  f 


(1  +  e  cos  uf 


sin(a;  +  u)  cos(cu  +  u)  sin  i 


A.  14.  Calculation  of  p  -/x 


df 

dM0 


The  Kepler  equation  (168)  entails  that 

df  1  df 

d  Ma  n  dt 

and,  therefore,  with  aid  of  (63)  we  obtain: 

df 


-.3/2 


-  /  X 


dM0 


VG 


,3/2 
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g 


—  w  JGma  (1  —  e2) 
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—  w  a 


2  VY 
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A. 15. 


The  terms 


/  _  r* 

\  d  (  _ 

(  M  X  f 

)  '  aa  G  *  f) 

In  the  course  of  numerical  computation,  it  is  convenient  to  keep  such  terms  as  parts  of 
the  Hamiltonian  (52).  In  case  these  terms  are  to  be  dealt  with  analytically,  one  will  need 
the  following  general  formulae  valid  for  any  Cj  ,  j  =  a  ,  e ,  12 ,  u  ,  i ,  Ma  : 


(185) 


51 


the  neglect  of  the  term  with  dft/d Cj  being  legitimate,  because  this  input  is  of  a  even 
higher  order  than  the  other  one.17  A  direct  algebraic  calculation  then  yields,  under  the  said 
approximation: 

(#  X  f)  '  ~QQ~  [p-  X  f)  =  (  l4  I  A{3  )  h  pjp  +  (  vt  +  )  f 2  +  (  hi  +  ^2  )  f 3 

(186) 

<9(7,  +  dC1,  J  ’ 


A*1  ^2 


,  dh  ,  dh  ,  \  ( 

hdc3  +  dcP2)  ~  ^  i 


/3  —  /^3  A6! 


expression  for  fj  given  above  by  (99  -  101) 

17  The  physical  meaning  of  dp,/dCj  is  the  influence  of  the  satellite  upon  the  precession  rate  of  the 
primary.  As  demonstrated  by  Laskar  (2004),  this  effect  is  of  a  higher  order  of  smallness. 
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